knitr::opts_chunk$set(message = FALSE, warning = FALSE)

require(data.table)
require(DESeq2)
if(!("qiime2R" %in% installed.packages()[,"Package"])){
    devtools::install_github("jbisanz/qiime2R")
}
require(qiime2R)
require(ggplot2)
require(vegan)
require(phyloseq)
require(umap)
require(MASS)
require(cowplot)
require(geosphere)
require(lme4)
require(nlme)
require(car)
require(jtools)

Statistical analysis

The main questions are:

* do microbiomes of Hydra from different locations differ
* do microbiomes differ between species
* do microbiomes differ between different sexual states 

The two common ways to test that is using an Anosim (a non parametric test on beta diversity) and an Adonis test (comparable to MANOVA on beta diversity). Common practice is to resample the libraries to the minimum number of reads per sample with the problem of loosing most of the data. However, the recommendation for statistical analysis based on the phyloseq authors are to normalize the reads with the DESeq2 package (which is also used to determine differential ASV between conditions) and the implementation of the variance stabilization transformation (vst). This requires, that at least one ASV is common to all samples, which is not the case for the diverse sample I got. To circumvent this problem, I will add a single read to every ASV in every sample. Another common and quite well working transformation is dividing by the library size. I will use the raw, the vst and the fraction to do the statistics.

otu <- read.csv(snakemake@input[["AbFilt_FT"]], row.names = 1) #"../data/dada2_AbFilt_FT.csv"
tax <- read.csv(snakemake@input[["AbFilt_Tax"]], row.names=1)#"../data/dada2_AbFilt_Tax.csv"
meta <- read.csv(snakemake@input[["AbFilt_meta"]], row.names=1)#"../data/dada2_AbFilt_meta.csv"
tree <- read_qza(snakemake@input[["AbFilt_fasta_tree"]])$data #"../data/dada2_AbFilt_fasta_tree.qza"

#otu <- read.csv("../data/dada2_AbFilt_FT.csv", row.names = 1)
#tax <- read.csv("../data/dada2_AbFilt_Tax.csv", row.names =1)
#meta <- read.csv("../data/dada2_AbFilt_meta.csv", row.names=1)
#tree <- read_qza("../data/dada2_AbFilt_fasta_tree.qza")$data

# remap the reproductive mode to the SEX, ASEX and NR
map_vect <- c("SEX","NR","SEX","ASEX","SEX","SEX","SEX","SEX")
names(map_vect) <- unique(meta$ReproductiveMode)
meta[["ReproductiveModeSimple"]] <- map_vect[meta[["ReproductiveMode"]]]


# artificially add a single count to all ESVs
otu2 <- otu+1

phy <- phyloseq(otu_table(otu, taxa_are_rows = TRUE),
                tax_table(as.matrix(tax)),
                sample_data(meta),
                phy_tree(tree))

otu_vst <- varianceStabilizingTransformation(as.matrix(otu2), fitType = "local")

# shift the vst correction to positive values
if(min(otu_vst) <0){
    otu_vst <- otu_vst + min(otu_vst)*-1
}

# create phyloseq objects
phy_vst <- phyloseq(otu_table(otu_vst,
                              taxa_are_rows = TRUE),
                    tax_table(as.matrix(tax)),
                    sample_data(meta),
                    phy_tree(tree))
phy_frac <- phyloseq(otu_table(otu/colSums(otu), taxa_are_rows = TRUE),
                tax_table(as.matrix(tax)),
                sample_data(meta),
                phy_tree(tree))

trnses <- list(raw = phy, vst = phy_vst, frac = phy_frac)

So lets look on the raw data here and whether we already see a separation there (on euclidean distance).

countPlots <- list()
plotDF <- data.frame()
for(n in names(trnses)){
    # calculate ordination
    pca <- prcomp(t(otu_table(trnses[[n]])))
    nmds <- isoMDS(dist(t(otu_table(trnses[[n]]))))
    ump <- umap(t(otu_table(trnses[[n]])))
    
    # create plots
    nPCA <- paste(n,"PCA",sep = "_")
    countPlots[[nPCA]] <- ggplot(as.data.frame(pca$x),
                                 aes(x=PC1, y=PC2, col = meta$PopID)) +
        geom_point() +
        ggtitle(nPCA) +
        theme(legend.position = "none")
    nNMDS <- paste(n,"NMDS", sep = "_")
    countPlots[[nNMDS]] <- ggplot(as.data.frame(nmds$points),
                                  aes(x=V1, y=V2, col = meta$PopID)) +
        geom_point() +
        ggtitle(nNMDS) +
        theme(legend.position = "none")
    nUMAP <- paste(n,"UMAP",sep ="_")
    countPlots[[nUMAP]] <- ggplot(as.data.frame(ump$layout),
                                  aes(x=V1, y=V2, col = meta$PopID))+
        geom_point() +
        ggtitle(nUMAP)
    
    
    # combine the plotting data in one data.frame
    plotDF <- rbind(plotDF,
                        cbind(V1 = pca$x[,1], V2 =pca$x[,2], ordination = "PCA", data = n, sample = rownames(pca$x)),
                        cbind(as.data.frame(nmds$points), ordination = "NMDS", data = n, sample = rownames(nmds$points)),
                        cbind(as.data.frame(ump$layout), ordination = "UMAP", data =n, sample = rownames(ump$layout)),
                    deparse.level = 2
    )
}
## initial  value 39.315036 
## final  value 39.315035 
## converged
## initial  value 37.542053 
## final  value 37.523047 
## converged
## initial  value 54.716981 
## iter   5 value 13.867358
## iter  10 value 12.441360
## iter  15 value 12.212251
## iter  20 value 12.078578
## iter  20 value 12.071296
## iter  20 value 12.071248
## final  value 12.071248 
## converged
# plot the different normalizations next to each other
plotDF<- merge(plotDF, meta, by.x="sample", by.y=0)
plotDF[,"V1"] <- as.numeric(plotDF[,"V1"])
plotDF[,"V2"] <- as.numeric(plotDF[,"V2"])

normPlots <- list()
for(n in names(trnses)){
    normPlots[[n]] <- ggplot(plotDF[plotDF[["data"]] == n,], aes(x = V1, y=V2, col =PopID)) +
            facet_wrap(.~ordination, scales = "free") +
            geom_point()+
            ggtitle(n) +
            theme(legend.position = "None")
}
cowplot::plot_grid(plotlist = normPlots, ncol =1)

To calculate differences between treatments, one needs to calculate \(\beta\)-diversity between samples. There are several distance measures, which can be calculated, the most common though are: Bray-Curtis, Jaccard, Unifrac and Weighted Unifrac. I will calculate all. To evaluate which might be the best to do statistics on, and which input data to use, I will plot ordination of the distance data. This is done by NMDS, PCoA, and UMAP.

## [1] "euclidean_raw"
## [1] "bray_raw"
## [1] "jaccard_raw"
## [1] "unifrac_raw"
## [1] "wunifrac_raw"
## [1] "euclidean_vst"
## [1] "bray_vst"
## [1] "jaccard_vst"
## [1] "unifrac_vst"
## [1] "wunifrac_vst"
## [1] "euclidean_frac"
## [1] "bray_frac"
## [1] "jaccard_frac"
## [1] "unifrac_frac"
## [1] "wunifrac_frac"
## [1] "euclidean_raw_NMDS"
## Run 0 stress 0.08369916 
## Run 1 stress 0.08493086 
## Run 2 stress 0.08730913 
## Run 3 stress 0.08644224 
## Run 4 stress 0.08512591 
## Run 5 stress 0.08413704 
## ... Procrustes: rmse 0.03578869  max resid 0.3719535 
## Run 6 stress 0.08379216 
## ... Procrustes: rmse 0.0224323  max resid 0.1555071 
## Run 7 stress 0.08985744 
## Run 8 stress 0.08678246 
## Run 9 stress 0.08613573 
## Run 10 stress 0.08426065 
## Run 11 stress 0.0838994 
## ... Procrustes: rmse 0.02361745  max resid 0.2184259 
## Run 12 stress 0.08653453 
## Run 13 stress 0.08513959 
## Run 14 stress 0.0867496 
## Run 15 stress 0.08454579 
## Run 16 stress 0.08578121 
## Run 17 stress 0.08490396 
## Run 18 stress 0.08705293 
## Run 19 stress 0.08369116 
## ... New best solution
## ... Procrustes: rmse 0.02715217  max resid 0.2166256 
## Run 20 stress 0.08508169 
## *** No convergence -- monoMDS stopping criteria:
##     18: no. of iterations >= maxit
##      2: stress ratio > sratmax
## [1] "euclidean_raw_PCoA"
## [1] "euclidean_raw_UMAP"
## [1] "bray_raw_NMDS"
## Run 0 stress 0.2773178 
## Run 1 stress 0.2811224 
## Run 2 stress 0.2814037 
## Run 3 stress 0.2797893 
## Run 4 stress 0.2781847 
## Run 5 stress 0.278967 
## Run 6 stress 0.2861727 
## Run 7 stress 0.2798618 
## Run 8 stress 0.2770567 
## ... New best solution
## ... Procrustes: rmse 0.02396036  max resid 0.1512107 
## Run 9 stress 0.2769995 
## ... New best solution
## ... Procrustes: rmse 0.02232634  max resid 0.1157627 
## Run 10 stress 0.2789943 
## Run 11 stress 0.2845254 
## Run 12 stress 0.2773256 
## ... Procrustes: rmse 0.0263711  max resid 0.1499273 
## Run 13 stress 0.3120813 
## Run 14 stress 0.2812445 
## Run 15 stress 0.2772958 
## ... Procrustes: rmse 0.02445731  max resid 0.1212653 
## Run 16 stress 0.2807581 
## Run 17 stress 0.277671 
## Run 18 stress 0.282163 
## Run 19 stress 0.2783355 
## Run 20 stress 0.2779285 
## *** No convergence -- monoMDS stopping criteria:
##      5: no. of iterations >= maxit
##     15: stress ratio > sratmax
## [1] "bray_raw_PCoA"
## [1] "bray_raw_UMAP"
## [1] "jaccard_raw_NMDS"
## Run 0 stress 0.2773784 
## Run 1 stress 0.2789441 
## Run 2 stress 0.2787378 
## Run 3 stress 0.3166571 
## Run 4 stress 0.2806936 
## Run 5 stress 0.280006 
## Run 6 stress 0.2770098 
## ... New best solution
## ... Procrustes: rmse 0.01917113  max resid 0.1242275 
## Run 7 stress 0.2842891 
## Run 8 stress 0.280299 
## Run 9 stress 0.2880282 
## Run 10 stress 0.2812404 
## Run 11 stress 0.2781231 
## Run 12 stress 0.2786286 
## Run 13 stress 0.2802902 
## Run 14 stress 0.2769373 
## ... New best solution
## ... Procrustes: rmse 0.02350488  max resid 0.155115 
## Run 15 stress 0.2806287 
## Run 16 stress 0.2821839 
## Run 17 stress 0.2811259 
## Run 18 stress 0.2792081 
## Run 19 stress 0.278262 
## Run 20 stress 0.2759648 
## ... New best solution
## ... Procrustes: rmse 0.01679744  max resid 0.1232613 
## *** No convergence -- monoMDS stopping criteria:
##      8: no. of iterations >= maxit
##     11: stress ratio > sratmax
##      1: scale factor of the gradient < sfgrmin
## [1] "jaccard_raw_PCoA"
## [1] "jaccard_raw_UMAP"
## [1] "unifrac_raw_NMDS"
## Run 0 stress 0.1950753 
## Run 1 stress 0.212001 
## Run 2 stress 0.200252 
## Run 3 stress 0.1996066 
## Run 4 stress 0.1976846 
## Run 5 stress 0.1995561 
## Run 6 stress 0.2033033 
## Run 7 stress 0.1950901 
## ... Procrustes: rmse 0.01680846  max resid 0.1909933 
## Run 8 stress 0.196432 
## Run 9 stress 0.1968841 
## Run 10 stress 0.199347 
## Run 11 stress 0.1949618 
## ... New best solution
## ... Procrustes: rmse 0.006644276  max resid 0.1063006 
## Run 12 stress 0.2136834 
## Run 13 stress 0.1962465 
## Run 14 stress 0.1961008 
## Run 15 stress 0.1969792 
## Run 16 stress 0.1991054 
## Run 17 stress 0.2014105 
## Run 18 stress 0.1950813 
## ... Procrustes: rmse 0.01574845  max resid 0.1912546 
## Run 19 stress 0.1949615 
## ... New best solution
## ... Procrustes: rmse 0.0005197987  max resid 0.006363654 
## ... Similar to previous best
## Run 20 stress 0.1961116 
## *** Solution reached
## [1] "unifrac_raw_PCoA"
## [1] "unifrac_raw_UMAP"
## [1] "wunifrac_raw_NMDS"
## Run 0 stress 0.161153 
## Run 1 stress 0.1818377 
## Run 2 stress 0.1772397 
## Run 3 stress 0.1875217 
## Run 4 stress 0.1737821 
## Run 5 stress 0.195064 
## Run 6 stress 0.1781208 
## Run 7 stress 0.1799434 
## Run 8 stress 0.1667898 
## Run 9 stress 0.176366 
## Run 10 stress 0.1945547 
## Run 11 stress 0.1856628 
## Run 12 stress 0.1836673 
## Run 13 stress 0.168293 
## Run 14 stress 0.1641538 
## Run 15 stress 0.1998411 
## Run 16 stress 0.1982581 
## Run 17 stress 0.1795297 
## Run 18 stress 0.1878017 
## Run 19 stress 0.1669929 
## Run 20 stress 0.1635855 
## *** No convergence -- monoMDS stopping criteria:
##     13: stress ratio > sratmax
##      7: scale factor of the gradient < sfgrmin
## [1] "wunifrac_raw_PCoA"
## [1] "wunifrac_raw_UMAP"
## [1] "euclidean_vst_NMDS"
## Run 0 stress 0.1609013 
## Run 1 stress 0.1626088 
## Run 2 stress 0.1648398 
## Run 3 stress 0.1649715 
## Run 4 stress 0.164929 
## Run 5 stress 0.1619001 
## Run 6 stress 0.1635389 
## Run 7 stress 0.1676298 
## Run 8 stress 0.1636399 
## Run 9 stress 0.1642341 
## Run 10 stress 0.165878 
## Run 11 stress 0.163316 
## Run 12 stress 0.1691651 
## Run 13 stress 0.1622519 
## Run 14 stress 0.1628481 
## Run 15 stress 0.1651573 
## Run 16 stress 0.1664558 
## Run 17 stress 0.1643859 
## Run 18 stress 0.1637193 
## Run 19 stress 0.1628497 
## Run 20 stress 0.1603392 
## ... New best solution
## ... Procrustes: rmse 0.02521003  max resid 0.1854854 
## *** No convergence -- monoMDS stopping criteria:
##     18: no. of iterations >= maxit
##      2: stress ratio > sratmax
## [1] "euclidean_vst_PCoA"
## [1] "euclidean_vst_UMAP"
## [1] "bray_vst_NMDS"
## Run 0 stress 0.2277878 
## Run 1 stress 0.2278856 
## ... Procrustes: rmse 0.007424575  max resid 0.0859814 
## Run 2 stress 0.2278728 
## ... Procrustes: rmse 0.009270762  max resid 0.0848336 
## Run 3 stress 0.2279974 
## ... Procrustes: rmse 0.01041228  max resid 0.1189395 
## Run 4 stress 0.2280305 
## ... Procrustes: rmse 0.00752595  max resid 0.08207802 
## Run 5 stress 0.2289447 
## Run 6 stress 0.2324784 
## Run 7 stress 0.2278391 
## ... Procrustes: rmse 0.0116298  max resid 0.08629948 
## Run 8 stress 0.2273921 
## ... New best solution
## ... Procrustes: rmse 0.01005542  max resid 0.08419402 
## Run 9 stress 0.2290124 
## Run 10 stress 0.2318222 
## Run 11 stress 0.2305865 
## Run 12 stress 0.2275014 
## ... Procrustes: rmse 0.009502425  max resid 0.0848474 
## Run 13 stress 0.2285419 
## Run 14 stress 0.2305028 
## Run 15 stress 0.417808 
## Run 16 stress 0.2295759 
## Run 17 stress 0.2280828 
## Run 18 stress 0.2303905 
## Run 19 stress 0.2288548 
## Run 20 stress 0.227982 
## *** No convergence -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     19: stress ratio > sratmax
## [1] "bray_vst_PCoA"
## [1] "bray_vst_UMAP"
## [1] "jaccard_vst_NMDS"
## Run 0 stress 0.227261 
## Run 1 stress 0.2328641 
## Run 2 stress 0.2273415 
## ... Procrustes: rmse 0.005575891  max resid 0.08523176 
## Run 3 stress 0.2290474 
## Run 4 stress 0.2273478 
## ... Procrustes: rmse 0.005451431  max resid 0.08526879 
## Run 5 stress 0.2305623 
## Run 6 stress 0.2317595 
## Run 7 stress 0.227494 
## ... Procrustes: rmse 0.003573452  max resid 0.0507231 
## Run 8 stress 0.2304353 
## Run 9 stress 0.2280415 
## Run 10 stress 0.2281836 
## Run 11 stress 0.2309814 
## Run 12 stress 0.2293575 
## Run 13 stress 0.2272952 
## ... Procrustes: rmse 0.006717131  max resid 0.083551 
## Run 14 stress 0.2324858 
## Run 15 stress 0.231942 
## Run 16 stress 0.2292307 
## Run 17 stress 0.2309612 
## Run 18 stress 0.2306358 
## Run 19 stress 0.2300717 
## Run 20 stress 0.227541 
## ... Procrustes: rmse 0.01020147  max resid 0.08705875 
## *** No convergence -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     19: stress ratio > sratmax
## [1] "jaccard_vst_PCoA"
## [1] "jaccard_vst_UMAP"
## [1] "unifrac_vst_NMDS"
## Run 0 stress 1.886457e-16 
## Run 1 stress 0 
## ... New best solution
## ... Procrustes: rmse 0.05733635  max resid 0.869531 
## Run 2 stress 0 
## ... Procrustes: rmse 0.06038091  max resid 0.2822453 
## Run 3 stress 0 
## ... Procrustes: rmse 0.06056925  max resid 0.2859596 
## Run 4 stress 0 
## ... Procrustes: rmse 0.06056342  max resid 0.2914502 
## Run 5 stress 0 
## ... Procrustes: rmse 0.06068206  max resid 0.3025266 
## Run 6 stress 0 
## ... Procrustes: rmse 0.06090321  max resid 0.3274585 
## Run 7 stress 0 
## ... Procrustes: rmse 0.06035597  max resid 0.2793788 
## Run 8 stress 0 
## ... Procrustes: rmse 0.06047073  max resid 0.300135 
## Run 9 stress 0 
## ... Procrustes: rmse 0.06069685  max resid 0.2949206 
## Run 10 stress 0 
## ... Procrustes: rmse 0.06040829  max resid 0.2830932 
## Run 11 stress 0 
## ... Procrustes: rmse 0.0601835  max resid 0.2831408 
## Run 12 stress 0 
## ... Procrustes: rmse 0.06082349  max resid 0.3174186 
## Run 13 stress 0 
## ... Procrustes: rmse 0.06002378  max resid 0.2948504 
## Run 14 stress 0 
## ... Procrustes: rmse 0.05956902  max resid 0.2573516 
## Run 15 stress 0 
## ... Procrustes: rmse 0.05996597  max resid 0.2811964 
## Run 16 stress 0 
## ... Procrustes: rmse 0.06054728  max resid 0.2931258 
## Run 17 stress 0 
## ... Procrustes: rmse 0.06052762  max resid 0.3052316 
## Run 18 stress 0 
## ... Procrustes: rmse 0.06018397  max resid 0.27517 
## Run 19 stress 0 
## ... Procrustes: rmse 0.06056692  max resid 0.2856627 
## Run 20 stress 0 
## ... Procrustes: rmse 0.06072305  max resid 0.3204826 
## *** No convergence -- monoMDS stopping criteria:
##     20: stress < smin
## [1] "unifrac_vst_PCoA"
## [1] "unifrac_vst_UMAP"
## [1] "wunifrac_vst_NMDS"
## Run 0 stress 0.1823021 
## Run 1 stress 0.417799 
## Run 2 stress 0.1823015 
## ... New best solution
## ... Procrustes: rmse 0.0005246753  max resid 0.007854011 
## ... Similar to previous best
## Run 3 stress 0.1858372 
## Run 4 stress 0.1823021 
## ... Procrustes: rmse 0.0004942914  max resid 0.007216553 
## ... Similar to previous best
## Run 5 stress 0.1858381 
## Run 6 stress 0.1866836 
## Run 7 stress 0.1861078 
## Run 8 stress 0.1823008 
## ... New best solution
## ... Procrustes: rmse 0.0002517982  max resid 0.003748713 
## ... Similar to previous best
## Run 9 stress 0.1865072 
## Run 10 stress 0.1823024 
## ... Procrustes: rmse 0.000312671  max resid 0.003092625 
## ... Similar to previous best
## Run 11 stress 0.1866844 
## Run 12 stress 0.4176636 
## Run 13 stress 0.1823022 
## ... Procrustes: rmse 0.0003406208  max resid 0.003818625 
## ... Similar to previous best
## Run 14 stress 0.1858385 
## Run 15 stress 0.1866838 
## Run 16 stress 0.1823019 
## ... Procrustes: rmse 0.0002926221  max resid 0.003797772 
## ... Similar to previous best
## Run 17 stress 0.1823026 
## ... Procrustes: rmse 0.0003126271  max resid 0.004146743 
## ... Similar to previous best
## Run 18 stress 0.1858369 
## Run 19 stress 0.1860942 
## Run 20 stress 0.1827007 
## ... Procrustes: rmse 0.004226909  max resid 0.06786374 
## *** Solution reached
## [1] "wunifrac_vst_PCoA"
## [1] "wunifrac_vst_UMAP"
## [1] "euclidean_frac_NMDS"
## Run 0 stress 0.1204247 
## Run 1 stress 0.1206017 
## ... Procrustes: rmse 0.0400446  max resid 0.3173344 
## Run 2 stress 0.1230261 
## Run 3 stress 0.1254943 
## Run 4 stress 0.1185844 
## ... New best solution
## ... Procrustes: rmse 0.04744261  max resid 0.3797871 
## Run 5 stress 0.121037 
## Run 6 stress 0.125874 
## Run 7 stress 0.1229305 
## Run 8 stress 0.1258877 
## Run 9 stress 0.1206976 
## Run 10 stress 0.1187527 
## ... Procrustes: rmse 0.04829711  max resid 0.4092316 
## Run 11 stress 0.1241809 
## Run 12 stress 0.1245341 
## Run 13 stress 0.1308419 
## Run 14 stress 0.1216227 
## Run 15 stress 0.1219218 
## Run 16 stress 0.1216191 
## Run 17 stress 0.1250719 
## Run 18 stress 0.1204038 
## Run 19 stress 0.1221799 
## Run 20 stress 0.1308524 
## *** No convergence -- monoMDS stopping criteria:
##     19: no. of iterations >= maxit
##      1: stress ratio > sratmax
## [1] "euclidean_frac_PCoA"
## [1] "euclidean_frac_UMAP"
## [1] "bray_frac_NMDS"
## Run 0 stress 0.2818453 
## Run 1 stress 0.2839507 
## Run 2 stress 0.2833467 
## Run 3 stress 0.2814406 
## ... New best solution
## ... Procrustes: rmse 0.01850348  max resid 0.1865741 
## Run 4 stress 0.2843496 
## Run 5 stress 0.2836991 
## Run 6 stress 0.2816992 
## ... Procrustes: rmse 0.0195696  max resid 0.1104788 
## Run 7 stress 0.2827069 
## Run 8 stress 0.2930529 
## Run 9 stress 0.283718 
## Run 10 stress 0.2852771 
## Run 11 stress 0.2856869 
## Run 12 stress 0.2837538 
## Run 13 stress 0.2819936 
## Run 14 stress 0.2857617 
## Run 15 stress 0.2827811 
## Run 16 stress 0.2871906 
## Run 17 stress 0.2831708 
## Run 18 stress 0.2845995 
## Run 19 stress 0.2855095 
## Run 20 stress 0.2811038 
## ... New best solution
## ... Procrustes: rmse 0.01760839  max resid 0.112224 
## *** No convergence -- monoMDS stopping criteria:
##      4: no. of iterations >= maxit
##     16: stress ratio > sratmax
## [1] "bray_frac_PCoA"
## [1] "bray_frac_UMAP"
## [1] "jaccard_frac_NMDS"
## Run 0 stress 0.2825359 
## Run 1 stress 0.295579 
## Run 2 stress 0.2864782 
## Run 3 stress 0.2829263 
## ... Procrustes: rmse 0.02673403  max resid 0.1903155 
## Run 4 stress 0.284254 
## Run 5 stress 0.2838242 
## Run 6 stress 0.2858169 
## Run 7 stress 0.2886027 
## Run 8 stress 0.2842911 
## Run 9 stress 0.2869574 
## Run 10 stress 0.2855451 
## Run 11 stress 0.2858612 
## Run 12 stress 0.2821634 
## ... New best solution
## ... Procrustes: rmse 0.02578425  max resid 0.1844668 
## Run 13 stress 0.2837911 
## Run 14 stress 0.2868005 
## Run 15 stress 0.2823502 
## ... Procrustes: rmse 0.02333661  max resid 0.1862586 
## Run 16 stress 0.282314 
## ... Procrustes: rmse 0.02790461  max resid 0.130169 
## Run 17 stress 0.2823868 
## ... Procrustes: rmse 0.02675134  max resid 0.1843313 
## Run 18 stress 0.2838576 
## Run 19 stress 0.28254 
## ... Procrustes: rmse 0.0246459  max resid 0.1156521 
## Run 20 stress 0.2843962 
## *** No convergence -- monoMDS stopping criteria:
##      7: no. of iterations >= maxit
##     13: stress ratio > sratmax
## [1] "jaccard_frac_PCoA"
## [1] "jaccard_frac_UMAP"
## [1] "unifrac_frac_NMDS"
## Run 0 stress 0.1950753 
## Run 1 stress 0.1963941 
## Run 2 stress 0.1960616 
## Run 3 stress 0.1969395 
## Run 4 stress 0.1982781 
## Run 5 stress 0.1955889 
## Run 6 stress 0.1969686 
## Run 7 stress 0.2170798 
## Run 8 stress 0.1947381 
## ... New best solution
## ... Procrustes: rmse 0.009979722  max resid 0.09647186 
## Run 9 stress 0.1977129 
## Run 10 stress 0.1988329 
## Run 11 stress 0.1984302 
## Run 12 stress 0.1996111 
## Run 13 stress 0.19627 
## Run 14 stress 0.1968442 
## Run 15 stress 0.1967855 
## Run 16 stress 0.197732 
## Run 17 stress 0.1999234 
## Run 18 stress 0.2007899 
## Run 19 stress 0.1968475 
## Run 20 stress 0.1978205 
## *** No convergence -- monoMDS stopping criteria:
##      1: no. of iterations >= maxit
##     16: stress ratio > sratmax
##      3: scale factor of the gradient < sfgrmin
## [1] "unifrac_frac_PCoA"
## [1] "unifrac_frac_UMAP"
## [1] "wunifrac_frac_NMDS"
## Run 0 stress 0.1693536 
## Run 1 stress 0.1822301 
## Run 2 stress 0.1862589 
## Run 3 stress 0.1724381 
## Run 4 stress 0.1700249 
## Run 5 stress 0.1971008 
## Run 6 stress 0.1890689 
## Run 7 stress 0.1929433 
## Run 8 stress 0.2079395 
## Run 9 stress 0.1728318 
## Run 10 stress 0.2044746 
## Run 11 stress 0.1960163 
## Run 12 stress 0.1776627 
## Run 13 stress 0.1936799 
## Run 14 stress 0.1741476 
## Run 15 stress 0.1897843 
## Run 16 stress 0.1836593 
## Run 17 stress 0.1833076 
## Run 18 stress 0.1710554 
## Run 19 stress 0.1721027 
## Run 20 stress 0.2050091 
## *** No convergence -- monoMDS stopping criteria:
##     12: stress ratio > sratmax
##      8: scale factor of the gradient < sfgrmin
## [1] "wunifrac_frac_PCoA"
## [1] "wunifrac_frac_UMAP"

From the plots it looks like the vst normalization in the UMAP clustering for euclidean distances gives best results for separation of Hydra from different locations directly. If one settles on one of the mostly used algorithms it looks like vst-bray and vst-jaccard give reliable results (at least UMAP is able to cluster them quite well). The data has many dimensions, which makes it difficult for NMDS an PCoA to reveal common patterns in the samples. Surprisingly the unifrac and wunifrac distances look rather scattered and give relatively bad clustering, which might imply that the inferred fast-tree ML-tree is no really representing the phylogenetic distance of the bacterial species, or that that the phylogenetic distance is not a good descriptor for the samples.

Lets see, whether the ordination is also representative for Species and Reproductive Mode differences.

sel <- paste(c("euclidean","bray","jaccard"),"vst", sep ="_")
plots2 <- list()
for(nm in sel){
    ordi <- umap(as.matrix(mtrxs[[nm]]), input="dist")
    plots2[[paste(nm,"RepMod",sep="_")]] <- ggplot(as.data.frame(ordi$layout), aes(x=V1,y=V2, color = meta$ReproductiveModeSimple))+
        ggtitle(nm) +
        labs(color = "RepMode")+
        geom_point()
    plots2[[paste(nm,"Species",sep="_")]] <- ggplot(as.data.frame(ordi$layout), aes(x=V1,y=V2, color = meta$Species))+
        ggtitle(nm) +
        labs(color = "Species")+
        geom_point()
    plots2[[paste(nm,"PopID",sep="_")]] <- ggplot(as.data.frame(ordi$layout), aes(x=V1,y=V2, color = meta$PopID))+
        ggtitle(nm) +
        labs(color = "PopID")+
        geom_point()
}
pp_UMAP <- cowplot::plot_grid(plotlist=plots2, ncol = 3)
pp_UMAP

Knowing these things, I will calculate an adonis on the euclidean, bray and jaccard distances for the vst normalized data and check for significant differences in sampling location (= Population ID [PopID]), reproductive mode (ReproductiveMode) and Species (Species).

sel <- paste(c("euclidean","bray","jaccard"),"vst", sep ="_")
stats <- list()
for(n in sel){
    stats[[n]] <- adonis(formula = mtrxs[[n]] ~meta$Species+meta$ReproductiveModeSimple+ meta$PopID, perm = 200)
}

print("Adonis testing for PopID")
## [1] "Adonis testing for PopID"
print(stats[sel])
## $euclidean_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ meta$Species + meta$ReproductiveModeSimple +      meta$PopID, permutations = 200) 
## 
## Permutation: free
## Number of permutations: 200
## 
## Terms added sequentially (first to last)
## 
##                              Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## meta$Species                  2     19538  9769.2  4.4897 0.02572 0.004975 **
## meta$ReproductiveModeSimple   2      8230  4115.1  1.8912 0.01083 0.004975 **
## meta$PopID                   20    209706 10485.3  4.8188 0.27604 0.004975 **
## Residuals                   240    522223  2175.9         0.68741            
## Total                       264    759698                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $bray_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ meta$Species + meta$ReproductiveModeSimple +      meta$PopID, permutations = 200) 
## 
## Permutation: free
## Number of permutations: 200
## 
## Terms added sequentially (first to last)
## 
##                              Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## meta$Species                  2     3.029 1.51455  7.9707 0.03920 0.004975 **
## meta$ReproductiveModeSimple   2     0.918 0.45886  2.4149 0.01188 0.004975 **
## meta$PopID                   20    27.732 1.38660  7.2974 0.35884 0.004975 **
## Residuals                   240    45.603 0.19001         0.59009            
## Total                       264    77.282                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $jaccard_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ meta$Species + meta$ReproductiveModeSimple +      meta$PopID, permutations = 200) 
## 
## Permutation: free
## Number of permutations: 200
## 
## Terms added sequentially (first to last)
## 
##                              Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## meta$Species                  2     3.019 1.50965  5.3374 0.03081 0.004975 **
## meta$ReproductiveModeSimple   2     1.030 0.51505  1.8210 0.01051 0.004975 **
## meta$PopID                   20    26.052 1.30260  4.6053 0.26588 0.004975 **
## Residuals                   240    67.883 0.28285         0.69279            
## Total                       264    97.984                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
stats0.5 <- list()
for(n in sel){
    stats0.5[[n]] <- adonis(formula = mtrxs[[n]] ~meta$PopID, perm = 200)
}
print("Adonis testing for PopID ONLY")
## [1] "Adonis testing for PopID ONLY"
print(stats0.5[sel])
## $euclidean_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ meta$PopID, permutations = 200) 
## 
## Permutation: free
## Number of permutations: 200
## 
## Terms added sequentially (first to last)
## 
##             Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## meta$PopID  20    223721 11186.1  5.0924 0.29449 0.004975 **
## Residuals  244    535977  2196.6         0.70551            
## Total      264    759698                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $bray_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ meta$PopID, permutations = 200) 
## 
## Permutation: free
## Number of permutations: 200
## 
## Terms added sequentially (first to last)
## 
##             Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## meta$PopID  20    29.950 1.49749  7.7196 0.38754 0.004975 **
## Residuals  244    47.332 0.19398         0.61246            
## Total      264    77.282                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $jaccard_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ meta$PopID, permutations = 200) 
## 
## Permutation: free
## Number of permutations: 200
## 
## Terms added sequentially (first to last)
## 
##             Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## meta$PopID  20    28.012 1.40060   4.884 0.28588 0.004975 **
## Residuals  244    69.972 0.28677         0.71412            
## Total      264    97.984                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
stats2 <- list()
for(n in sel){
    stats2[[n]] <- adonis(formula = mtrxs[[n]] ~  meta$Species, strata= meta$PopID, perm =200)
}

print("Adonis testing for Species and Reproductive Mode")
## [1] "Adonis testing for Species and Reproductive Mode"
print(stats2)
## $euclidean_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ meta$Species, permutations = 200,      strata = meta$PopID) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 200
## 
## Terms added sequentially (first to last)
## 
##               Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## meta$Species   2     19538  9769.2  3.4581 0.02572 0.004975 **
## Residuals    262    740160  2825.0         0.97428            
## Total        264    759698                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $bray_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ meta$Species, permutations = 200,      strata = meta$PopID) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 200
## 
## Terms added sequentially (first to last)
## 
##               Df SumsOfSqs MeanSqs F.Model     R2   Pr(>F)   
## meta$Species   2     3.029 1.51455  5.3441 0.0392 0.004975 **
## Residuals    262    74.253 0.28341         0.9608            
## Total        264    77.282                 1.0000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $jaccard_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ meta$Species, permutations = 200,      strata = meta$PopID) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 200
## 
## Terms added sequentially (first to last)
## 
##               Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## meta$Species   2     3.019 1.50965   4.165 0.03081 0.004975 **
## Residuals    262    94.965 0.36246         0.96919            
## Total        264    97.984                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# select only oligactis for species test in adonis

sel2 <- rownames(meta[meta$Species == "OLI",])

stats3 <- list()
for(n in sel){
    stats3[[n]] <- adonis(formula = as.matrix(mtrxs[[n]])[sel2,sel2] ~  meta[sel2,"ReproductiveModeSimple"], strata= meta[sel2,"PopID"], perm =200)
}
print("Adonis testing for Reproductive Mode")
## [1] "Adonis testing for Reproductive Mode"
print(stats3)
## $euclidean_vst
## 
## Call:
## adonis(formula = as.matrix(mtrxs[[n]])[sel2, sel2] ~ meta[sel2,      "ReproductiveModeSimple"], permutations = 200, strata = meta[sel2,      "PopID"]) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 200
## 
## Terms added sequentially (first to last)
## 
##                                       Df SumsOfSqs MeanSqs F.Model      R2
## meta[sel2, "ReproductiveModeSimple"]   2      8536  4268.0   1.376 0.01318
## Residuals                            206    638971  3101.8         0.98682
## Total                                208    647507                 1.00000
##                                      Pr(>F)
## meta[sel2, "ReproductiveModeSimple"] 0.4826
## Residuals                                  
## Total                                      
## 
## $bray_vst
## 
## Call:
## adonis(formula = as.matrix(mtrxs[[n]])[sel2, sel2] ~ meta[sel2,      "ReproductiveModeSimple"], permutations = 200, strata = meta[sel2,      "PopID"]) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 200
## 
## Terms added sequentially (first to last)
## 
##                                       Df SumsOfSqs MeanSqs F.Model      R2
## meta[sel2, "ReproductiveModeSimple"]   2     0.973 0.48635  1.6038 0.01533
## Residuals                            206    62.468 0.30324         0.98467
## Total                                208    63.441                 1.00000
##                                      Pr(>F)
## meta[sel2, "ReproductiveModeSimple"]  0.408
## Residuals                                  
## Total                                      
## 
## $jaccard_vst
## 
## Call:
## adonis(formula = as.matrix(mtrxs[[n]])[sel2, sel2] ~ meta[sel2,      "ReproductiveModeSimple"], permutations = 200, strata = meta[sel2,      "PopID"]) 
## 
## Blocks:  strata 
## Permutation: free
## Number of permutations: 200
## 
## Terms added sequentially (first to last)
## 
##                                       Df SumsOfSqs MeanSqs F.Model      R2
## meta[sel2, "ReproductiveModeSimple"]   2     1.062 0.53086  1.4001 0.01341
## Residuals                            206    78.109 0.37917         0.98659
## Total                                208    79.171                 1.00000
##                                      Pr(>F)
## meta[sel2, "ReproductiveModeSimple"] 0.4478
## Residuals                                  
## Total

To finally answer the main question from the beginning: 1.) There is a consistent and rather strong signature in the data which is indicative for the sampling location, hence bacterial communities in Hydra differ if they live in different locations. 2.) There is a consistent difference in the bacterial communities from different Hydra species (this is what was expected). 3.) There is a less strong, but detectable change in the microbial community if the animals have different reproductive modes. This is only detectable if the model design is stratified by the sampling location of the 16S data. However, it indicates, that change in the reproductive mode is associated with the differences in the microbiota.

We obtained new meta data for the lakes where the animals were sampled and Sebastian suggested to test for differences in the microbiota in dependence of the water body. Does it make a difference if Hydra lives in a lake or a flowing (river) water body? Further we obtained the geocoordinates of the sampling sites. I will use the distance between geolocation and \(\beta\)-diversity to correlate these two frameworks.

metaLakes <- read.csv(snakemake@input[["MetaLakes"]])

# join the meta data
meta2 <- merge(meta, metaLakes, by.x = "PopID", by.y = "SiteID", all.x=TRUE)


stats3 <- list()
for(n in sel){
    stats3[[n]] <- adonis(formula = mtrxs[[n]] ~ meta2$waterbody, perm =500)
}
stats3
## $euclidean_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ meta2$waterbody, permutations = 500) 
## 
## Permutation: free
## Number of permutations: 500
## 
## Terms added sequentially (first to last)
## 
##                  Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## meta2$waterbody   1      9838  9838.3  3.4506 0.01295 0.001996 **
## Residuals       263    749860  2851.2         0.98705            
## Total           264    759698                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $bray_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ meta2$waterbody, permutations = 500) 
## 
## Permutation: free
## Number of permutations: 500
## 
## Terms added sequentially (first to last)
## 
##                  Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## meta2$waterbody   1     1.110 1.10971  3.8315 0.01436 0.001996 **
## Residuals       263    76.172 0.28963         0.98564            
## Total           264    77.282                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $jaccard_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ meta2$waterbody, permutations = 500) 
## 
## Permutation: free
## Number of permutations: 500
## 
## Terms added sequentially (first to last)
## 
##                  Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## meta2$waterbody   1     1.082 1.08223  2.9373 0.01104 0.001996 **
## Residuals       263    96.902 0.36845         0.98896            
## Total           264    97.984                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# The waterbody is a factor which changes the microbiota!

stats4 <- list()
for(n in sel){
    stats4[[n]] <- adonis(formula = mtrxs[[n]] ~ factor(meta2$NutrientLoad), perm =500)
}
stats4
## $euclidean_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ factor(meta2$NutrientLoad), permutations = 500) 
## 
## Permutation: free
## Number of permutations: 500
## 
## Terms added sequentially (first to last)
## 
##                             Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## factor(meta2$NutrientLoad)   2     12479  6239.4  2.1877 0.01643 0.001996 **
## Residuals                  262    747219  2852.0         0.98357            
## Total                      264    759698                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $bray_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ factor(meta2$NutrientLoad), permutations = 500) 
## 
## Permutation: free
## Number of permutations: 500
## 
## Terms added sequentially (first to last)
## 
##                             Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## factor(meta2$NutrientLoad)   2     1.694 0.84709  2.9362 0.02192 0.001996 **
## Residuals                  262    75.588 0.28850         0.97808            
## Total                      264    77.282                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## $jaccard_vst
## 
## Call:
## adonis(formula = mtrxs[[n]] ~ factor(meta2$NutrientLoad), permutations = 500) 
## 
## Permutation: free
## Number of permutations: 500
## 
## Terms added sequentially (first to last)
## 
##                             Df SumsOfSqs MeanSqs F.Model      R2   Pr(>F)   
## factor(meta2$NutrientLoad)   2     1.697 0.84834  2.3083 0.01732 0.001996 **
## Residuals                  262    96.288 0.36751         0.98268            
## Total                      264    97.984                 1.00000            
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# calculate the geographic distance between sampling points
geoDist <- distm(metaLakes[,c("lon","lat")])
colnames(geoDist) <- metaLakes[["SiteID"]]
rownames(geoDist) <- metaLakes[["SiteID"]]

# set the upper triangle to NA, so its easily filtered afterwards
geoDist[upper.tri(geoDist, diag=FALSE)] <- NA

# create a data frame which associates distance of microbiome to geographic distance
geoDistLong <- data.table(reshape2::melt(geoDist))
geoDistLong <- geoDistLong[!is.na(value),]

mtrxsLong <- list()
for(mat in sel){
    bacDist <- as.matrix(mtrxs[[mat]])
    bacDist[upper.tri(bacDist, diag=TRUE)] <- NA
    bacDistLong <- data.table(reshape2::melt(bacDist))
    bacDistLong <- bacDistLong[!is.na(value),]
    # add the PopID information to associate it back to the geographic distance
    bacDistLong <- cbind(bacDistLong, 
                          Var1PopID = meta[bacDistLong[,Var1], "PopID"],
                          Var2PopID = meta[bacDistLong[,Var2], "PopID"],
                          GeoDist = -1)
    mtrxsLong[[mat]] <- bacDistLong
}

# merge both tables
for(i in 1:nrow(geoDistLong)){
    for(mat in names(mtrxsLong)){
       sel2 <-  ((mtrxsLong[[mat]][,Var1PopID] %in% geoDistLong[i,Var1] &
                 mtrxsLong[[mat]][,Var2PopID] %in% geoDistLong[i,Var2]) |
                (mtrxsLong[[mat]][,Var2PopID] %in% geoDistLong[i,Var1] &
                 mtrxsLong[[mat]][,Var1PopID] %in% geoDistLong[i,Var2]))
        mtrxsLong[[mat]][sel2, "GeoDist"] <- geoDistLong[i,value]
    }
}


# create simple dotplots
geoDistPlots <- list()
geoDistModels <- list()
for(mat in names(mtrxsLong)){
    nm <- strsplit(mat, split="_")[[1]]
    p <- ggplot(mtrxsLong[[mat]], aes(y=value, x=GeoDist)) +
        geom_point(alpha = 0.1) +
        xlab("Geographic distance") +
        ylab(paste(nm[1], "distance")) +
        ggtitle(paste("Distance based on", nm[2]))
    geoDistPlots[[mat]] <- p
    print(cor.test(method ="pearson", mtrxsLong[[mat]][,value], mtrxsLong[[mat]][,GeoDist]))
    geoDistModels[[mat]] <- lm(data=mtrxsLong[[mat]], value~GeoDist)
}
## 
##  Pearson's product-moment correlation
## 
## data:  mtrxsLong[[mat]][, value] and mtrxsLong[[mat]][, GeoDist]
## t = -14.565, df = 34978, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.08804867 -0.06721599
## sample estimates:
##         cor 
## -0.07764081 
## 
## 
##  Pearson's product-moment correlation
## 
## data:  mtrxsLong[[mat]][, value] and mtrxsLong[[mat]][, GeoDist]
## t = 25.377, df = 34978, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1241526 0.1447327
## sample estimates:
##       cor 
## 0.1344572 
## 
## 
##  Pearson's product-moment correlation
## 
## data:  mtrxsLong[[mat]][, value] and mtrxsLong[[mat]][, GeoDist]
## t = 27.216, df = 34978, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.1337289 0.1542533
## sample estimates:
##       cor 
## 0.1440066
pp <- cowplot::plot_grid(plotlist = geoDistPlots, ncol = 1)
pp

# create a distance matrix for the distance in geo-location of the same dimensions like the beta-diversity to perform a mantel test
# 1. cast into wide format
geoDist_Samples <- dcast(mtrxsLong[[1]], Var1~Var2, Value.var = "GeoDist")
# 2. format as data.frame, correct rownames, which are missing in data.table
geoDist_Samples_names <- geoDist_Samples[, Var1]
geoDist_Samples <- as.data.frame(geoDist_Samples)
rownames(geoDist_Samples) <- geoDist_Samples_names
geoDist_Samples <- geoDist_Samples[,-1]
# 3. add an additional column and row, because I left out the diagonal by the first conversion from distance to long format
geoDist_Samples <- cbind(geoDist_Samples, added = NA)
colnames(geoDist_Samples)[ncol(geoDist_Samples)] <- geoDist_Samples_names[length(geoDist_Samples_names)]
geoDist_Samples_names <- colnames(geoDist_Samples)
geoDist_Samples <- rbind(added = NA, geoDist_Samples)
rownames(geoDist_Samples)[1] <- geoDist_Samples_names[1]
# 4. reformat into distance object
geoDist_Samples <- as.dist(geoDist_Samples)


# perform a mantel test on the distance matrices
for(mat in names(mtrxsLong)){
    print(mat)
    print(mantel(geoDist_Samples, mtrxs[[mat]]))
}
## [1] "euclidean_vst"
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = geoDist_Samples, ydis = mtrxs[[mat]]) 
## 
## Mantel statistic r:     1 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0770 0.0981 0.1168 0.1451 
## Permutation: free
## Number of permutations: 999
## 
## [1] "bray_vst"
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = geoDist_Samples, ydis = mtrxs[[mat]]) 
## 
## Mantel statistic r: 0.7914 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0553 0.0727 0.0867 0.1004 
## Permutation: free
## Number of permutations: 999
## 
## [1] "jaccard_vst"
## 
## Mantel statistic based on Pearson's product-moment correlation 
## 
## Call:
## mantel(xdis = geoDist_Samples, ydis = mtrxs[[mat]]) 
## 
## Mantel statistic r: 0.7649 
##       Significance: 0.001 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0583 0.0726 0.0860 0.0959 
## Permutation: free
## Number of permutations: 999

So there is a high statistic significant association of geographic distance and \(\beta\)-diversity, however, the effect is marginal and probably only significant because of the shear amount of data in the set. Its a nice feature, but I would not make to much a deal out of it.

# The idea of Jacint was to associate the alpha diversity to the 
alphaDiv <- estimate_richness(phy)
alphaNames <- colnames(alphaDiv)
alphaDiv <- merge(alphaDiv, meta, by=0, all.x=TRUE)
# transform the data to fit normality a little better
alphaDiv[["Simpson"]] <- alphaDiv[["Simpson"]]^2
alphaDiv[["Chao1"]] <- log(alphaDiv[["Chao1"]])


alphaPlots <- list()
for(nm in c("Simpson","Shannon","Chao1")){
    for(group in c("PopID", "Species", "ReproductiveModeSimple")){
        id <- paste(nm, group, sep ="_")
        p <- ggplot(alphaDiv, aes_string(x = group, y= nm, fill = group))+
                    geom_boxplot() +
                    ylab(nm) +
                    xlab(group) +
                    theme(legend.position = "None")
        alphaPlots[[id]] <- p
    }
}
pp <- cowplot::plot_grid(plotlist = alphaPlots, rel_widths = rep(c(3,1,1),3))

print(pp)

There seems to be major differences in the alpha diversity of the different populations, but also in the different species. Lets fit some linear (mixed) models to the data and see if these differences are significant.

Lets start with the Simpson index:

# lets fit a linear mixed model to the data and check for population
simpFull <- lm(data=alphaDiv, Simpson ~ PopID+Species+ReproductiveModeSimple)
simpPopSpec <- lm(data=alphaDiv, Simpson ~ PopID+Species)
simpPopNull <- lm(data=alphaDiv, Simpson ~ PopID)
simpPopLMM <- lme(data=alphaDiv, Simpson ~ PopID, random = ~1|Species/ReproductiveModeSimple, method = "ML")

anova(simpPopLMM, simpFull)
##            Model df       AIC      BIC   logLik   Test  L.Ratio p-value
## simpPopLMM     1 24 -22.26486 63.64866 35.13243                        
## simpFull       2 26 -32.50277 60.57021 42.25138 1 vs 2 14.23791   8e-04
anova(simpPopLMM, simpPopNull)
##             Model df       AIC      BIC   logLik   Test  L.Ratio p-value
## simpPopLMM      1 24 -22.26486 63.64866 35.13243                        
## simpPopNull     2 22 -15.98438 62.76967 29.99219 1 vs 2 10.28047  0.0059
# the full model is still better, test the difference to the pop+species model
anova(simpPopNull,simpPopSpec,simpFull)
## Analysis of Variance Table
## 
## Model 1: Simpson ~ PopID
## Model 2: Simpson ~ PopID + Species
## Model 3: Simpson ~ PopID + Species + ReproductiveModeSimple
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    244 12.373                                  
## 2    242 11.553  2   0.81973 8.7210 0.0002207 ***
## 3    240 11.279  2   0.27367 2.9115 0.0563182 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(simpPopNull,simpPopSpec,simpFull)
##             df       AIC
## simpPopNull 22 -15.98438
## simpPopSpec 24 -30.14997
## simpFull    26 -32.50277
BIC(simpPopNull,simpPopSpec,simpFull)
##             df      BIC
## simpPopNull 22 62.76967
## simpPopSpec 24 55.76355
## simpFull    26 60.57021
# the linear model with only Population and Species is best fitting the data

# lets check for species differences
simpSpecNull <- lm(data=alphaDiv, Simpson ~ Species)
simpSpecLMM1 <- lme(data=alphaDiv, Simpson ~ Species, random =~1|PopID, method = "ML")
simpSpecLMM2 <- lme(data=alphaDiv, Simpson ~ Species+ReproductiveModeSimple, random =~1|PopID, method = "ML")

anova(simpSpecLMM2, simpSpecLMM1)
##              Model df       AIC      BIC   logLik   Test  L.Ratio p-value
## simpSpecLMM2     1  7 -5.466983 19.59113 9.733491                        
## simpSpecLMM1     2  5 -4.822704 13.07595 7.411352 1 vs 2 4.644279  0.0981
anova(simpSpecLMM2, simpFull)
##              Model df       AIC      BIC   logLik   Test  L.Ratio p-value
## simpSpecLMM2     1  7  -5.46698 19.59113  9.73349                        
## simpFull         2 26 -32.50277 60.57021 42.25138 1 vs 2 65.03578  <.0001
anova(simpSpecLMM2, simpSpecNull)
##              Model df      AIC      BIC     logLik   Test  L.Ratio p-value
## simpSpecLMM2     1  7 -5.46698 19.59113   9.733491                        
## simpSpecNull     2  4 36.93405 51.25297 -14.467024 1 vs 2 48.40103  <.0001
anova(simpSpecNull, simpPopSpec, simpFull)
## Analysis of Variance Table
## 
## Model 1: Simpson ~ Species
## Model 2: Simpson ~ PopID + Species
## Model 3: Simpson ~ PopID + Species + ReproductiveModeSimple
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    262 17.306                                  
## 2    242 11.553 20    5.7527 6.1203 4.551e-13 ***
## 3    240 11.279  2    0.2737 2.9115   0.05632 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(simpSpecNull,simpPopSpec,simpFull)
##              df       AIC
## simpSpecNull  4  36.93405
## simpPopSpec  24 -30.14997
## simpFull     26 -32.50277
BIC(simpSpecNull,simpPopSpec,simpFull)
##              df      BIC
## simpSpecNull  4 51.25297
## simpPopSpec  24 55.76355
## simpFull     26 60.57021
# again the population+Species simple linear model is best fitting the data

simpMod <- simpPopSpec

# ok lets repeat this procedure with the other diversity measures
shanFull <- lm(data=alphaDiv, Shannon ~ PopID+Species+ReproductiveModeSimple)
shanPopSpec <- lm(data=alphaDiv, Shannon ~ PopID+Species)
shanPopNull <- lm(data= alphaDiv, Shannon ~ PopID)
shanSpecNull <- lm(data= alphaDiv, Shannon ~ Species)
shanPopLMM <- lme(data=alphaDiv, Shannon ~ PopID, random = ~1|Species/ReproductiveModeSimple, method = "ML")
shanSpecLMM1 <- lme(data=alphaDiv, Shannon~ Species, random =~1|PopID, method = "ML")
shanSpecLMM2 <- lme(data=alphaDiv, Shannon~ Species+ReproductiveModeSimple, random =~1|PopID, method = "ML")

# find the best population model
anova(shanPopLMM, shanPopNull)
##             Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## shanPopLMM      1 24 634.9987 720.9122 -293.4993                        
## shanPopNull     2 22 638.2704 717.0244 -297.1352 1 vs 2 7.271693  0.0264
anova(shanPopLMM, shanFull)
##            Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## shanPopLMM     1 24 634.9987 720.9122 -293.4993                        
## shanFull       2 26 625.7518 718.8248 -286.8759 1 vs 2 13.24687  0.0013
anova(shanPopNull, shanPopSpec, shanFull)
## Analysis of Variance Table
## 
## Model 1: Shannon ~ PopID
## Model 2: Shannon ~ PopID + Species
## Model 3: Shannon ~ PopID + Species + ReproductiveModeSimple
##   Res.Df    RSS Df Sum of Sq      F   Pr(>F)   
## 1    244 146.11                                
## 2    242 139.26  2    6.8489 6.0777 0.002662 **
## 3    240 135.23  2    4.0375 3.5829 0.029292 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(shanPopNull, shanPopSpec, shanFull)
##             df      AIC
## shanPopNull 22 638.2704
## shanPopSpec 24 629.5482
## shanFull    26 625.7518
BIC(shanPopNull, shanPopSpec, shanFull)
##             df      BIC
## shanPopNull 22 717.0244
## shanPopSpec 24 715.4617
## shanFull    26 718.8248
# this one is a bit tricky, but I would go to the full model here, as log-likelihood and AIC are in concordance here


# species effect
anova(shanSpecLMM1, shanSpecLMM2)
##              Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## shanSpecLMM1     1  5 655.3448 673.2434 -322.6724                        
## shanSpecLMM2     2  7 653.1819 678.2400 -319.5910 1 vs 2 6.162849  0.0459
anova(shanSpecLMM2, shanSpecNull)
##              Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## shanSpecLMM2     1  7 653.1819 678.2400 -319.5910                        
## shanSpecNull     2  4 707.2619 721.5808 -349.6309 1 vs 2 60.07994  <.0001
anova(shanSpecLMM2, shanFull)
##              Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## shanSpecLMM2     1  7 653.1819 678.2400 -319.5910                        
## shanFull         2 26 625.7518 718.8248 -286.8759 1 vs 2 65.43015  <.0001
anova(shanSpecNull, shanPopSpec, shanFull)
## Analysis of Variance Table
## 
## Model 1: Shannon ~ Species
## Model 2: Shannon ~ PopID + Species
## Model 3: Shannon ~ PopID + Species + ReproductiveModeSimple
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    262 217.15                                  
## 2    242 139.26 20    77.883 6.9113 5.385e-15 ***
## 3    240 135.23  2     4.038 3.5829   0.02929 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(shanSpecNull, shanPopSpec, shanFull)
##              df      AIC
## shanSpecNull  4 707.2619
## shanPopSpec  24 629.5482
## shanFull     26 625.7518
BIC(shanSpecNull, shanPopSpec, shanFull)
##              df      BIC
## shanSpecNull  4 721.5808
## shanPopSpec  24 715.4617
## shanFull     26 718.8248
# again the full model seems to be the best, at least concord in log-likelihood and AIC

shanMod <- shanFull


chaoFull <- lm(data=alphaDiv, Chao1 ~ PopID+Species+ReproductiveModeSimple)
chaoPopNull <- lm(data=alphaDiv, Chao1 ~ PopID)
chaoPopSpec <- lm(data=alphaDiv, Chao1 ~ PopID+Species)
chaoSpecNull<- lm(data=alphaDiv, Chao1 ~ Species)
chaoPopLMM <- lme(data=alphaDiv, Chao1 ~ PopID, random = ~1|Species/ReproductiveModeSimple, method = "ML")
chaoSpecLMM1 <- lme(data=alphaDiv, Chao1 ~ Species, random =~1|PopID, method = "ML")
chaoSpecLMM2 <- lme(data=alphaDiv, Chao1 ~ Species+ReproductiveModeSimple, random =~1|PopID, method = "ML")

# population
anova(chaoPopLMM, chaoPopNull)
##             Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## chaoPopLMM      1 24 447.1783 533.0919 -199.5892                        
## chaoPopNull     2 22 451.5532 530.3073 -203.7766 1 vs 2 8.374895  0.0152
anova(chaoPopLMM, chaoFull)
##            Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## chaoPopLMM     1 24 447.1783 533.0919 -199.5892                        
## chaoFull       2 26 438.2494 531.3224 -193.1247 1 vs 2 12.92893  0.0016
anova(chaoPopNull, chaoPopSpec, chaoFull)
## Analysis of Variance Table
## 
## Model 1: Chao1 ~ PopID
## Model 2: Chao1 ~ PopID + Species
## Model 3: Chao1 ~ PopID + Species + ReproductiveModeSimple
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    244 72.225                                  
## 2    242 67.124  2    5.1012 9.1851 0.0001433 ***
## 3    240 66.646  2    0.4778 0.8603 0.4243201    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(chaoPopNull, chaoPopSpec, chaoFull)
##             df      AIC
## chaoPopNull 22 451.5532
## chaoPopSpec 24 436.1425
## chaoFull    26 438.2494
BIC(chaoPopNull, chaoPopSpec, chaoFull)
##             df      BIC
## chaoPopNull 22 530.3073
## chaoPopSpec 24 522.0560
## chaoFull    26 531.3224
# the linear model with population and species is fitting the data best
chaoMod <- chaoPopSpec

# species
anova(chaoSpecLMM1,chaoSpecLMM2 , chaoSpecNull)
##              Model df      AIC      BIC    logLik   Test   L.Ratio p-value
## chaoSpecLMM1     1  5 477.4263 495.3250 -233.7132                         
## chaoSpecLMM2     2  7 479.7240 504.7821 -232.8620 1 vs 2   1.70232  0.4269
## chaoSpecNull     3  4 580.3695 594.6884 -286.1848 2 vs 3 106.64548  <.0001
anova(chaoSpecLMM1,chaoSpecLMM2 , chaoFull)
##              Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## chaoSpecLMM1     1  5 477.4263 495.3250 -233.7132                        
## chaoSpecLMM2     2  7 479.7240 504.7821 -232.8620 1 vs 2  1.70232  0.4269
## chaoFull         3 26 438.2494 531.3224 -193.1247 2 vs 3 79.47462  <.0001
anova(chaoSpecNull, chaoPopSpec, chaoFull)
## Analysis of Variance Table
## 
## Model 1: Chao1 ~ Species
## Model 2: Chao1 ~ PopID + Species
## Model 3: Chao1 ~ PopID + Species + ReproductiveModeSimple
##   Res.Df     RSS Df Sum of Sq       F Pr(>F)    
## 1    262 134.523                                
## 2    242  67.124 20    67.399 12.1356 <2e-16 ***
## 3    240  66.646  2     0.478  0.8603 0.4243    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC(chaoSpecNull, chaoPopSpec, chaoFull)
##              df      AIC
## chaoSpecNull  4 580.3695
## chaoPopSpec  24 436.1425
## chaoFull     26 438.2494
BIC(chaoSpecNull, chaoPopSpec, chaoFull)
##              df      BIC
## chaoSpecNull  4 594.6884
## chaoPopSpec  24 522.0560
## chaoFull     26 531.3224

OK after fitting the data to a (rather simple) model lets see if the assumptions hold an plot some diagnostics.

simp_pp <- plot_grid(~plot(simpMod, which =1),
                     ~plot(simpMod, which =2),
                     ~plot(simpMod,which=4),
                     ~plot(simpMod,which=5))

shan_pp <- plot_grid(~plot(shanMod, which =1),
                     ~plot(shanMod, which =2),
                     ~plot(shanMod,which=4),
                     ~plot(shanMod,which=5))

chao_pp <- plot_grid(~plot(chaoMod, which =1),
                     ~plot(chaoMod, which =2),
                     ~plot(chaoMod,which=4),
                     ~plot(chaoMod,which=5))
simp_pp

shan_pp

chao_pp

The models look quite okay, lets look on the coefficients.

summary(simpMod)
## 
## Call:
## lm(formula = Simpson ~ PopID + Species, data = alphaDiv)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.63243 -0.12389  0.01434  0.13754  0.51192 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.18747    0.09809   1.911 0.057170 .  
## PopIDM108    0.02174    0.09559   0.227 0.820279    
## PopIDM109    0.27082    0.09783   2.768 0.006072 ** 
## PopIDM26     0.18345    0.10050   1.825 0.069189 .  
## PopIDM28     0.20177    0.08539   2.363 0.018920 *  
## PopIDM31     0.17151    0.10769   1.593 0.112574    
## PopIDM34     0.19550    0.11293   1.731 0.084699 .  
## PopIDM44     0.19147    0.08821   2.171 0.030928 *  
## PopIDM67    -0.36262    0.09771  -3.711 0.000256 ***
## PopIDM70    -0.08486    0.10815  -0.785 0.433406    
## PopIDM71    -0.09660    0.09203  -1.050 0.294922    
## PopIDM72     0.03747    0.08408   0.446 0.656265    
## PopIDM78    -0.18167    0.16931  -1.073 0.284339    
## PopIDM79     0.07531    0.08931   0.843 0.399930    
## PopIDM83    -0.03842    0.09090  -0.423 0.672861    
## PopIDM84     0.20800    0.09203   2.260 0.024696 *  
## PopIDM85     0.26310    0.08933   2.945 0.003541 ** 
## PopIDM86     0.55867    0.15980   3.496 0.000562 ***
## PopIDM89     0.13784    0.08809   1.565 0.118969    
## PopIDM90     0.16183    0.08287   1.953 0.051991 .  
## PopIDR12    -0.07562    0.09783  -0.773 0.440283    
## SpeciesOLI   0.26641    0.07092   3.756 0.000216 ***
## SpeciesVUL   0.29167    0.07246   4.025 7.62e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2185 on 242 degrees of freedom
## Multiple R-squared:  0.3467, Adjusted R-squared:  0.2873 
## F-statistic: 5.836 on 22 and 242 DF,  p-value: 3.653e-13
Anova(simpMod)
## Anova Table (Type II tests)
## 
## Response: Simpson
##            Sum Sq  Df F value    Pr(>F)    
## PopID      5.7527  20  6.0251 7.434e-13 ***
## Species    0.8197   2  8.5854 0.0002499 ***
## Residuals 11.5530 242                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(shanMod)
## 
## Call:
## lm(formula = Shannon ~ PopID + Species + ReproductiveModeSimple, 
##     data = alphaDiv)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.46721 -0.47374  0.00618  0.40648  1.92327 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 0.4848     0.3661   1.324 0.186773    
## PopIDM108                   0.1324     0.3345   0.396 0.692450    
## PopIDM109                   1.1456     0.3361   3.408 0.000766 ***
## PopIDM26                    1.1531     0.3530   3.266 0.001250 ** 
## PopIDM28                    0.7970     0.3009   2.649 0.008607 ** 
## PopIDM31                    0.6990     0.3738   1.870 0.062734 .  
## PopIDM34                    0.5082     0.3883   1.309 0.191765    
## PopIDM44                    0.7452     0.3109   2.397 0.017309 *  
## PopIDM67                   -1.1852     0.3445  -3.441 0.000684 ***
## PopIDM70                   -0.1457     0.3733  -0.390 0.696635    
## PopIDM71                   -0.3625     0.3246  -1.117 0.265133    
## PopIDM72                    0.4837     0.2958   1.635 0.103288    
## PopIDM78                   -0.6623     0.5864  -1.129 0.259867    
## PopIDM79                    0.2035     0.3114   0.654 0.514057    
## PopIDM83                    0.1703     0.3223   0.528 0.597673    
## PopIDM84                    0.6658     0.3256   2.045 0.041942 *  
## PopIDM85                    1.4690     0.3124   4.703 4.33e-06 ***
## PopIDM86                    1.4996     0.5514   2.720 0.007011 ** 
## PopIDM89                    0.6154     0.3090   1.991 0.047575 *  
## PopIDM90                    0.4680     0.2883   1.623 0.105886    
## PopIDR12                   -0.2334     0.3363  -0.694 0.488338    
## SpeciesOLI                  1.1089     0.2634   4.210 3.62e-05 ***
## SpeciesVUL                  1.0162     0.2750   3.695 0.000272 ***
## ReproductiveModeSimpleNR    0.1121     0.1493   0.751 0.453500    
## ReproductiveModeSimpleSEX   0.3651     0.1364   2.677 0.007944 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7506 on 240 degrees of freedom
## Multiple R-squared:  0.3936, Adjusted R-squared:  0.3329 
## F-statistic:  6.49 on 24 and 240 DF,  p-value: 1.109e-15
Anova(shanMod)
## Anova Table (Type II tests)
## 
## Response: Shannon
##                         Sum Sq  Df F value    Pr(>F)    
## PopID                   80.756  20  7.1663 1.319e-15 ***
## Species                 10.113   2  8.9744 0.0001743 ***
## ReproductiveModeSimple   4.038   2  3.5829 0.0292918 *  
## Residuals              135.227 240                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(chaoMod)
## 
## Call:
## lm(formula = Chao1 ~ PopID + Species, data = alphaDiv)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.71014 -0.30002  0.02297  0.30643  1.05792 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.28383    0.23645  13.888  < 2e-16 ***
## PopIDM108    0.44896    0.23040   1.949 0.052496 .  
## PopIDM109    0.61489    0.23581   2.608 0.009686 ** 
## PopIDM26     0.72523    0.24226   2.994 0.003042 ** 
## PopIDM28     0.48852    0.20582   2.373 0.018404 *  
## PopIDM31     0.05747    0.25959   0.221 0.824972    
## PopIDM34    -0.42703    0.27221  -1.569 0.118010    
## PopIDM44     0.28525    0.21261   1.342 0.180971    
## PopIDM67    -0.86959    0.23553  -3.692 0.000275 ***
## PopIDM70    -0.75293    0.26069  -2.888 0.004225 ** 
## PopIDM71    -0.83685    0.22182  -3.773 0.000203 ***
## PopIDM72     0.44506    0.20267   2.196 0.029045 *  
## PopIDM78    -1.79411    0.40811  -4.396 1.65e-05 ***
## PopIDM79     0.18240    0.21528   0.847 0.397694    
## PopIDM83     0.15393    0.21909   0.703 0.482985    
## PopIDM84     0.42829    0.22182   1.931 0.054679 .  
## PopIDM85     0.94649    0.21531   4.396 1.65e-05 ***
## PopIDM86     0.55613    0.38518   1.444 0.150081    
## PopIDM89     0.36098    0.21234   1.700 0.090415 .  
## PopIDM90     0.01105    0.19975   0.055 0.955916    
## PopIDR12    -0.57231    0.23581  -2.427 0.015955 *  
## SpeciesOLI   0.68218    0.17095   3.991 8.74e-05 ***
## SpeciesVUL   0.37232    0.17466   2.132 0.034044 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5267 on 242 degrees of freedom
## Multiple R-squared:  0.5293, Adjusted R-squared:  0.4865 
## F-statistic: 12.37 on 22 and 242 DF,  p-value: < 2.2e-16
Anova(chaoMod)
## Anova Table (Type II tests)
## 
## Response: Chao1
##           Sum Sq  Df F value    Pr(>F)    
## PopID     67.399  20 12.1496 < 2.2e-16 ***
## Species    5.101   2  9.1957 0.0001415 ***
## Residuals 67.124 242                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There is a clear effect of species, and population on alpha diversity in the samples. The populations M78 and M67 have the lowest alpha diversity, while M85 shows highest diversity in the models and plots. In the Shannon index, there is also some statistic difference for the reproductive mode. Sexually propagating animals show the highest alpha diversity. This coincides with the expression of periculin in the egg fleck and might be a consequence of the destroyed homoeostasis of the microbiota in these animals.

Finally we can have a look on the pairwise comparisons.

Anova(simpMod)
## Anova Table (Type II tests)
## 
## Response: Simpson
##            Sum Sq  Df F value    Pr(>F)    
## PopID      5.7527  20  6.0251 7.434e-13 ***
## Species    0.8197   2  8.5854 0.0002499 ***
## Residuals 11.5530 242                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(shanMod)
## Anova Table (Type II tests)
## 
## Response: Shannon
##                         Sum Sq  Df F value    Pr(>F)    
## PopID                   80.756  20  7.1663 1.319e-15 ***
## Species                 10.113   2  8.9744 0.0001743 ***
## ReproductiveModeSimple   4.038   2  3.5829 0.0292918 *  
## Residuals              135.227 240                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(chaoMod)
## Anova Table (Type II tests)
## 
## Response: Chao1
##           Sum Sq  Df F value    Pr(>F)    
## PopID     67.399  20 12.1496 < 2.2e-16 ***
## Species    5.101   2  9.1957 0.0001415 ***
## Residuals 67.124 242                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
TukeyHSD(aov(simpMod))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = simpMod)
## 
## $PopID
##                   diff          lwr          upr     p adj
## M108-M107  0.019213394 -0.325400158  0.363826945 1.0000000
## M109-M107  0.268291064 -0.084432152  0.621014279 0.4220948
## M26-M107   0.180923214 -0.181465441  0.543311868 0.9671343
## M28-M107   0.210073212 -0.092960012  0.513106436 0.6074779
## M31-M107   0.172587803 -0.216094250  0.561269855 0.9909131
## M34-M107   0.192973661 -0.214316025  0.600263348 0.9810422
## M44-M107   0.188941744 -0.128998666  0.506882153 0.8478043
## M67-M107  -0.362616483 -0.715339698 -0.009893268 0.0361681
## M70-M107  -0.068445582 -0.442565048  0.305673884 1.0000000
## M71-M107  -0.099122430 -0.430872703  0.232627843 0.9999643
## M72-M107   0.034943603 -0.268089621  0.337976827 1.0000000
## M78-M107  -0.184197789 -0.795132319  0.426736740 0.9999588
## M79-M107   0.056710004 -0.265280765  0.378700773 1.0000000
## M83-M107  -0.033733892 -0.360292135  0.292824350 1.0000000
## M84-M107   0.205475656 -0.126274617  0.537225928 0.7960575
## M85-M107   0.260570771 -0.061419998  0.582561541 0.3029158
## M86-M107   0.289732328 -0.229462186  0.808926843 0.9075285
## M89-M107   0.042561710 -0.250921745  0.336045165 1.0000000
## M90-M107   0.121734394 -0.173375021  0.416843808 0.9963106
## R12-M107  -0.078146873 -0.430870088  0.274576342 0.9999998
## M109-M108  0.249077670 -0.095535882  0.593691221 0.5251314
## M26-M108   0.161709820 -0.192790349  0.516209990 0.9875813
## M28-M108   0.190859818 -0.102694211  0.484413847 0.7226320
## M31-M108   0.153374409 -0.227963461  0.534712279 0.9973342
## M34-M108   0.173760268 -0.226526782  0.574047317 0.9930626
## M44-M108   0.169728350 -0.139190632  0.478647332 0.9190472
## M67-M108  -0.381829876 -0.726443428 -0.037216325 0.0133830
## M70-M108  -0.087658976 -0.454142547  0.278824596 0.9999991
## M71-M108  -0.118335824 -0.441450473  0.204778826 0.9992486
## M72-M108   0.015730209 -0.277823821  0.309284238 1.0000000
## M78-M108  -0.203411183 -0.809699756  0.402877389 0.9997883
## M79-M108   0.037496610 -0.275589462  0.350582683 1.0000000
## M83-M108  -0.052947286 -0.370728845  0.264834272 1.0000000
## M84-M108   0.186262262 -0.136852388  0.509376912 0.8794462
## M85-M108   0.241357378 -0.071728695  0.554443450 0.3953962
## M86-M108   0.270518935 -0.243200617  0.784238486 0.9447680
## M89-M108   0.023348316 -0.260337023  0.307033655 1.0000000
## M90-M108   0.102521000 -0.182846132  0.387888132 0.9994273
## R12-M108  -0.097360267 -0.441973819  0.247253284 0.9999855
## M26-M109  -0.087367850 -0.449756504  0.275020805 0.9999990
## M28-M109  -0.058217852 -0.361251076  0.244815372 1.0000000
## M31-M109  -0.095703261 -0.484385313  0.292978792 0.9999985
## M34-M109  -0.075317402 -0.482607089  0.331972284 1.0000000
## M44-M109  -0.079349320 -0.397289729  0.238591090 0.9999981
## M67-M109  -0.630907546 -0.983630761 -0.278184331 0.0000001
## M70-M109  -0.336736645 -0.710856111  0.037382821 0.1409852
## M71-M109  -0.367413494 -0.699163766 -0.035663221 0.0134696
## M72-M109  -0.233347461 -0.536380685  0.069685763 0.3976013
## M78-M109  -0.452488853 -1.063423383  0.158445677 0.4757642
## M79-M109  -0.211581060 -0.533571829  0.110409710 0.7048733
## M83-M109  -0.302024956 -0.628583198  0.024533286 0.1109854
## M84-M109  -0.062815408 -0.394565681  0.268934865 1.0000000
## M85-M109  -0.007720292 -0.329711061  0.314270477 1.0000000
## M86-M109   0.021441265 -0.497753250  0.540635779 1.0000000
## M89-M109  -0.225729354 -0.519212809  0.067754101 0.3999205
## M90-M109  -0.146556670 -0.441666084  0.148552745 0.9688563
## R12-M109  -0.346437937 -0.699161152  0.006285278 0.0609831
## M28-M26    0.029149998 -0.285080817  0.343380813 1.0000000
## M31-M26   -0.008335411 -0.405809446  0.389138624 1.0000000
## M34-M26    0.012050448 -0.403637848  0.427738743 1.0000000
## M44-M26    0.008018530 -0.320611923  0.336648983 1.0000000
## M67-M26   -0.543539697 -0.905928351 -0.181151042 0.0000298
## M70-M26   -0.249368796 -0.632614468  0.133876876 0.7213473
## M71-M26   -0.280045644 -0.622054582  0.061963294 0.2818894
## M72-M26   -0.145979611 -0.460210426  0.168251204 0.9847417
## M78-M26   -0.365121003 -0.981686385  0.251444378 0.8517284
## M79-M26   -0.124213210 -0.456763846  0.208337427 0.9990118
## M83-M26   -0.214657106 -0.551632114  0.122317902 0.7551789
## M84-M26    0.024552442 -0.317456496  0.366561380 1.0000000
## M85-M26    0.079647558 -0.252903079  0.412198194 0.9999991
## M86-M26    0.108809115 -0.416999610  0.634617839 0.9999999
## M89-M26   -0.138361504 -0.443393320  0.166670312 0.9883579
## M90-M26   -0.059188820 -0.365785358  0.247407717 1.0000000
## R12-M26   -0.259070087 -0.621458742  0.103318567 0.5470988
## M31-M28   -0.037485409 -0.381708020  0.306737202 1.0000000
## M34-M28   -0.017099550 -0.382202764  0.348003664 1.0000000
## M44-M28   -0.021131468 -0.282859518  0.240596582 1.0000000
## M67-M28   -0.572689695 -0.875722918 -0.269656471 0.0000000
## M70-M28   -0.278518794 -0.606208957  0.049171370 0.2199190
## M71-M28   -0.309195642 -0.587536641 -0.030854643 0.0128919
## M72-M28   -0.175129609 -0.418531752  0.068272534 0.5342313
## M78-M28   -0.394271001 -0.977928836  0.189386833 0.6560621
## M79-M28   -0.153363208 -0.419996896  0.113270480 0.8815440
## M83-M28   -0.243807104 -0.515938973  0.028324764 0.1465833
## M84-M28   -0.004597556 -0.282938555  0.273743443 1.0000000
## M85-M28    0.050497560 -0.216136129  0.317131248 1.0000000
## M86-M28    0.079659117 -0.407145169  0.566463402 1.0000000
## M89-M28   -0.167511502 -0.398915924  0.063892920 0.5220242
## M90-M28   -0.088338818 -0.321801952  0.145124316 0.9988192
## R12-M28   -0.288220085 -0.591253309  0.014813139 0.0850220
## M34-M31    0.020385859 -0.418413595  0.459185312 1.0000000
## M44-M31    0.016353941 -0.341062029  0.373769911 1.0000000
## M67-M31   -0.535204286 -0.923886338 -0.146522233 0.0002452
## M70-M31   -0.241033385 -0.649231188  0.167164418 0.8548950
## M71-M31   -0.271710233 -0.641464631  0.098044165 0.4917075
## M72-M31   -0.137644200 -0.481866812  0.206578411 0.9975287
## M78-M31   -0.356785592 -0.989162909  0.275591725 0.8986697
## M79-M31   -0.115877799 -0.476901519  0.245145922 0.9998913
## M83-M31   -0.206321695 -0.571424909  0.158781519 0.8973185
## M84-M31    0.032887853 -0.336866545  0.402642251 1.0000000
## M85-M31    0.087982969 -0.273040752  0.449006689 0.9999988
## M86-M31    0.117144526 -0.427119212  0.661408263 0.9999999
## M89-M31   -0.130026093 -0.465872203  0.205820017 0.9983889
## M90-M31   -0.050853409 -0.388121312  0.286414493 1.0000000
## R12-M31   -0.250734676 -0.639416729  0.137947376 0.7353870
## M44-M34   -0.004031918 -0.381599562  0.373535727 1.0000000
## M67-M34   -0.555590144 -0.962879831 -0.148300458 0.0003028
## M70-M34   -0.261419243 -0.687372993  0.164534506 0.8083787
## M71-M34   -0.292096092 -0.681363950  0.097171767 0.4493933
## M72-M34   -0.158030059 -0.523133273  0.207073155 0.9933006
## M78-M34   -0.377171451 -1.021152989  0.266810087 0.8635955
## M79-M34   -0.136263657 -0.517248273  0.244720959 0.9994625
## M83-M34   -0.226707554 -0.611560133  0.158145025 0.8574778
## M84-M34    0.012501994 -0.376765864  0.401769853 1.0000000
## M85-M34    0.067597110 -0.313387506  0.448581726 1.0000000
## M86-M34    0.096758667 -0.460945705  0.654463039 1.0000000
## M89-M34   -0.150411952 -0.507628638  0.206804734 0.9952237
## M90-M34   -0.071239268 -0.429793014  0.287314478 1.0000000
## R12-M34   -0.271120535 -0.678410221  0.136169152 0.6824325
## M67-M44   -0.551558226 -0.869498636 -0.233617817 0.0000004
## M70-M44   -0.257387326 -0.598910110  0.084135459 0.4406173
## M71-M44   -0.288064174 -0.582564906  0.006436558 0.0637383
## M72-M44   -0.153998141 -0.415726192  0.107729909 0.8587294
## M78-M44   -0.373139533 -0.964674348  0.218395281 0.7699602
## M79-M44   -0.132231740 -0.415693365  0.151229885 0.9840219
## M83-M44   -0.222675636 -0.511315071  0.065963798 0.3939389
## M84-M44    0.016533912 -0.277966820  0.311034644 1.0000000
## M85-M44    0.071629028 -0.211832598  0.355090653 0.9999977
## M86-M44    0.100790585 -0.395430517  0.597011686 0.9999999
## M89-M44   -0.146380034 -0.396989244  0.104229175 0.8664106
## M90-M44   -0.067207350 -0.319718744  0.185304044 0.9999945
## R12-M44   -0.267088617 -0.585029027  0.050851792 0.2385305
## M70-M67    0.294170901 -0.079948565  0.668290367 0.3567695
## M71-M67    0.263494053 -0.068256220  0.595244325 0.3375773
## M72-M67    0.397560085  0.094526861  0.700593309 0.0007001
## M78-M67    0.178418693 -0.432515836  0.789353223 0.9999752
## M79-M67    0.419326487  0.097335718  0.741317256 0.0008146
## M83-M67    0.328882590  0.002324348  0.655440833 0.0461147
## M84-M67    0.568092138  0.236341866  0.899842411 0.0000006
## M85-M67    0.623187254  0.301196485  0.945178023 0.0000000
## M86-M67    0.652348811  0.133154296  1.171543326 0.0016478
## M89-M67    0.405178192  0.111694737  0.698661647 0.0002310
## M90-M67    0.484350876  0.189241462  0.779460291 0.0000022
## R12-M67    0.284469609 -0.068253606  0.637192825 0.3091257
## M71-M70   -0.030676848 -0.385091791  0.323738095 1.0000000
## M72-M70    0.103389184 -0.224300979  0.431079348 0.9999165
## M78-M70   -0.115752208 -0.739284651  0.507780236 1.0000000
## M79-M70    0.125155586 -0.220141046  0.470452217 0.9993506
## M83-M70    0.034711690 -0.314848025  0.384271404 1.0000000
## M84-M70    0.273921238 -0.080493705  0.628336181 0.3903185
## M85-M70    0.329016353 -0.016280278  0.674312985 0.0834989
## M86-M70    0.358177910 -0.175783373  0.892139194 0.6687808
## M89-M70    0.111007291 -0.207872383  0.429886966 0.9996358
## M90-M70    0.190179975 -0.130196795  0.510556746 0.8490435
## R12-M70   -0.009701291 -0.383820757  0.364418175 1.0000000
## M72-M71    0.134066033 -0.144274966  0.412407032 0.9773652
## M78-M71   -0.085075359 -0.684145953  0.513995234 1.0000000
## M79-M71    0.155832434 -0.143036488  0.454701356 0.9496754
## M83-M71    0.065388538 -0.238395699  0.369172775 0.9999999
## M84-M71    0.304598086 -0.004760639  0.613956810 0.0593734
## M85-M71    0.359693201  0.060824280  0.658562123 0.0035989
## M86-M71    0.388854758 -0.116325923  0.894035440 0.3983862
## M89-M71    0.141684140 -0.126228375  0.409596654 0.9425275
## M90-M71    0.220856824 -0.048835860  0.490549507 0.2816857
## R12-M71    0.020975557 -0.310774716  0.352725829 1.0000000
## M78-M72   -0.219141392 -0.802799227  0.364516442 0.9989384
## M79-M72    0.021766401 -0.244867287  0.288400090 1.0000000
## M83-M72   -0.068677495 -0.340809363  0.203454374 0.9999978
## M84-M72    0.170532053 -0.107808946  0.448873052 0.8106371
## M85-M72    0.225627169 -0.041006519  0.492260857 0.2268384
## M86-M72    0.254788726 -0.232015560  0.741593011 0.9478359
## M89-M72    0.007618107 -0.223786315  0.239022529 1.0000000
## M90-M72    0.086790791 -0.146672343  0.320253925 0.9990744
## R12-M72   -0.113090476 -0.416123700  0.189942748 0.9990235
## M79-M78    0.240907794 -0.352813849  0.834629436 0.9970124
## M83-M78    0.150463897 -0.445747155  0.746674949 0.9999978
## M84-M78    0.389673445 -0.209397149  0.988744039 0.7218870
## M85-M78    0.444768561 -0.148953081  1.038490203 0.4527683
## M86-M78    0.473930118 -0.246063130  1.193923366 0.7019185
## M89-M78    0.226759499 -0.351997679  0.805516677 0.9981107
## M90-M78    0.305932183 -0.273651201  0.885515568 0.9435248
## R12-M78    0.106050916 -0.504883613  0.716985446 1.0000000
## M83-M79   -0.090443896 -0.383538889  0.202651096 0.9999408
## M84-M79    0.148765652 -0.150103270  0.447634573 0.9681169
## M85-M79    0.203860767 -0.084136532  0.491858067 0.5667594
## M86-M79    0.233022324 -0.265803630  0.731848279 0.9837722
## M89-M79   -0.014148294 -0.269876524  0.241579935 1.0000000
## M90-M79    0.065024390 -0.192568226  0.322617005 0.9999978
## R12-M79   -0.134856877 -0.456847647  0.187133892 0.9955318
## M84-M83    0.239209548 -0.064574689  0.542993785 0.3540089
## M85-M83    0.294304664  0.001209671  0.587399656 0.0477141
## M86-M83    0.323466221 -0.178320151  0.825252593 0.7365308
## M89-M83    0.076295602 -0.185160240  0.337751444 0.9999755
## M90-M83    0.155468286 -0.107811386  0.418747958 0.8548499
## R12-M83   -0.044412981 -0.370971223  0.282145261 1.0000000
## M85-M84    0.055095116 -0.243773806  0.353964037 1.0000000
## M86-M84    0.084256673 -0.420924009  0.589437354 1.0000000
## M89-M84   -0.162913946 -0.430826460  0.104998568 0.8203072
## M90-M84   -0.083741262 -0.353433946  0.185951421 0.9999348
## R12-M84   -0.283622529 -0.615372802  0.048127744 0.2109243
## M86-M85    0.029161557 -0.469664398  0.527987512 1.0000000
## M89-M85   -0.218009062 -0.473737291  0.037719168 0.2152699
## M90-M85   -0.138836378 -0.396428993  0.118756238 0.9316658
## R12-M85   -0.338717645 -0.660708414 -0.016726876 0.0271592
## M89-M86   -0.247170619 -0.728088300  0.233747062 0.9561666
## M90-M86   -0.167997935 -0.649909591  0.313913721 0.9996283
## R12-M86   -0.367879202 -0.887073716  0.151315313 0.5648014
## M90-M89    0.079172684 -0.141753635  0.300099003 0.9994472
## R12-M89   -0.120708583 -0.414192038  0.172774872 0.9964458
## R12-M90   -0.199881267 -0.494990681  0.095228148 0.6511985
## 
## $Species
##               diff         lwr       upr     p adj
## OLI-CIR 0.16630613  0.03973661 0.2928756 0.0061370
## VUL-CIR 0.20972749  0.06229557 0.3571594 0.0026458
## VUL-OLI 0.04342137 -0.04744624 0.1342890 0.4984479
TukeyHSD(aov(shanMod))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = shanMod)
## 
## $PopID
##                   diff          lwr         upr     p adj
## M108-M107  0.268667819 -0.915352007  1.45268764 0.9999996
## M109-M107  1.154821697 -0.057061231  2.36670462 0.0834317
## M26-M107   1.232022859 -0.013068489  2.47711421 0.0562408
## M28-M107   0.736880899 -0.304277649  1.77803945 0.5668822
## M31-M107   0.722758564 -0.612671381  2.05818851 0.9289949
## M34-M107   0.481000927 -0.918360942  1.88036280 0.9996965
## M44-M107   0.775958473 -0.316418035  1.86833498 0.5596522
## M67-M107  -1.047537142 -2.259420070  0.16434579 0.1944955
## M70-M107  -0.256454693 -1.541850648  1.02894126 1.0000000
## M71-M107  -0.417760280 -1.557584580  0.72206402 0.9992386
## M72-M107   0.527413943 -0.513744605  1.56857249 0.9619185
## M78-M107  -0.633481997 -2.732524801  1.46556081 0.9999581
## M79-M107   0.264682481 -0.841610214  1.37097517 0.9999991
## M83-M107   0.086633488 -1.035352085  1.20861906 1.0000000
## M84-M107   0.582482032 -0.557342268  1.72230633 0.9584811
## M85-M107   1.400606061  0.294313366  2.50689875 0.0014275
## M86-M107   0.672486082 -1.111357287  2.45632945 0.9988747
## M89-M107   0.219014597 -0.789332950  1.22736214 0.9999998
## M90-M107   0.410004549 -0.603929453  1.42393855 0.9971351
## R12-M107  -0.235372697 -1.447255625  0.97651023 1.0000000
## M109-M108  0.886153878 -0.297865947  2.07017370 0.4544826
## M26-M108   0.963355040 -0.254633130  2.18134321 0.3453938
## M28-M108   0.468213081 -0.540376945  1.47680311 0.9848373
## M31-M108   0.454090745 -0.856106131  1.76428762 0.9996575
## M34-M108   0.212333109 -1.162969170  1.58763539 1.0000000
## M44-M108   0.507290655 -0.554090122  1.56867143 0.9791278
## M67-M108  -1.316204961 -2.500224786 -0.13218514 0.0127654
## M70-M108  -0.525122512 -1.784283136  0.73403811 0.9957549
## M71-M108  -0.686428099 -1.796582210  0.42372601 0.7982246
## M72-M108   0.258746124 -0.749843901  1.26733615 0.9999970
## M78-M108  -0.902149816 -2.985230085  1.18093045 0.9932401
## M79-M108  -0.003985338 -1.079683364  1.07171269 1.0000000
## M83-M108  -0.182034331 -1.273865061  0.90979640 1.0000000
## M84-M108   0.313814214 -0.796339898  1.42396833 0.9999853
## M85-M108   1.131938242  0.056240216  2.20763627 0.0270541
## M86-M108   0.403818263 -1.361214281  2.16885081 0.9999996
## M89-M108  -0.049653222 -1.024336499  0.92503006 1.0000000
## M90-M108   0.141336731 -0.839124833  1.12179829 1.0000000
## R12-M108  -0.504040516 -1.688060341  0.67997931 0.9945188
## M26-M109   0.077201162 -1.167890186  1.32229251 1.0000000
## M28-M109  -0.417940797 -1.459099346  0.62321775 0.9973948
## M31-M109  -0.432063133 -1.767493078  0.90336681 0.9998770
## M34-M109  -0.673820770 -2.073182639  0.72554110 0.9773982
## M44-M109  -0.378863223 -1.471239732  0.71351329 0.9996540
## M67-M109  -2.202358839 -3.414241767 -0.99047591 0.0000001
## M70-M109  -1.411276390 -2.696672345 -0.12588044 0.0152573
## M71-M109  -1.572581977 -2.712406277 -0.43275768 0.0002351
## M72-M109  -0.627407754 -1.668566302  0.41375079 0.8313677
## M78-M109  -1.788303694 -3.887346498  0.31073911 0.2161813
## M79-M109  -0.890139216 -1.996431911  0.21615348 0.3133016
## M83-M109  -1.068188209 -2.190173782  0.05379736 0.0841867
## M84-M109  -0.572339665 -1.712163965  0.56748464 0.9651186
## M85-M109   0.245784364 -0.860508331  1.35207706 0.9999998
## M86-M109  -0.482335615 -2.266178984  1.30150775 0.9999929
## M89-M109  -0.935807100 -1.944154647  0.07254045 0.1075074
## M90-M109  -0.744817148 -1.758751150  0.26911685 0.4922593
## R12-M109  -1.390194394 -2.602077322 -0.17831147 0.0079417
## M28-M26   -0.495141960 -1.574773081  0.58448916 0.9867739
## M31-M26   -0.509264295 -1.874901645  0.85637306 0.9990311
## M34-M26   -0.751021932 -2.179239659  0.67719580 0.9454193
## M44-M26   -0.456064386 -1.585169639  0.67304087 0.9971764
## M67-M26   -2.279560001 -3.524651349 -1.03446865 0.0000001
## M70-M26   -1.488477553 -2.805229235 -0.17172587 0.0099540
## M71-M26   -1.649783139 -2.824854060 -0.47471222 0.0001568
## M72-M26   -0.704608916 -1.784240037  0.37502220 0.7161787
## M78-M26   -1.865504857 -3.983894086  0.25288437 0.1685547
## M79-M26   -0.967340379 -2.109914560  0.17523380 0.2259716
## M83-M26   -1.145389371 -2.303164765  0.01238602 0.0563672
## M84-M26   -0.649540827 -1.824611748  0.52553009 0.9146032
## M85-M26    0.168583201 -0.973990980  1.31115738 1.0000000
## M86-M26   -0.559536777 -2.366105182  1.24703163 0.9999371
## M89-M26   -1.013008262 -2.061033553  0.03501703 0.0722077
## M90-M26   -0.822018310 -1.875419656  0.23138304 0.3713590
## R12-M26   -1.467395556 -2.712486904 -0.22230421 0.0051405
## M31-M28   -0.014122335 -1.196798973  1.16855430 1.0000000
## M34-M28   -0.255879972 -1.510297978  0.99853803 0.9999999
## M44-M28    0.039077574 -0.860165068  0.93832022 1.0000000
## M67-M28   -1.784418041 -2.825576590 -0.74325949 0.0000005
## M70-M28   -0.993335593 -2.119210209  0.13253902 0.1660609
## M71-M28   -1.154641180 -2.110962426 -0.19831993 0.0034074
## M72-M28   -0.209466957 -1.045745627  0.62681171 0.9999980
## M78-M28   -1.370362897 -3.375688702  0.63496291 0.6345224
## M79-M28   -0.472198419 -1.388295803  0.44389897 0.9548575
## M83-M28   -0.650247412 -1.585235389  0.28474057 0.6011272
## M84-M28   -0.154398867 -1.110720114  0.80192238 1.0000000
## M85-M28    0.663725161 -0.252372223  1.57982255 0.5201264
## M86-M28   -0.064394818 -1.736952158  1.60816252 1.0000000
## M89-M28   -0.517866302 -1.312923324  0.27719072 0.7194437
## M90-M28   -0.326876350 -1.129006672  0.47525397 0.9968349
## R12-M28   -0.972253596 -2.013412145  0.06890495 0.1014822
## M34-M31   -0.241757637 -1.749380450  1.26586518 1.0000000
## M44-M31    0.053199909 -1.174806338  1.28120616 1.0000000
## M67-M31   -1.770295706 -3.105725651 -0.43486576 0.0005647
## M70-M31   -0.979213258 -2.381695224  0.42326871 0.5935126
## M71-M31   -1.140518844 -2.410917341  0.12987965 0.1440394
## M72-M31   -0.195344621 -1.378021259  0.98733202 1.0000000
## M78-M31   -1.356240562 -3.528956281  0.81647516 0.7850136
## M79-M31   -0.458076084 -1.698477806  0.78232564 0.9991537
## M83-M31   -0.636125076 -1.890543082  0.61829293 0.9615119
## M84-M31   -0.140276532 -1.410675029  1.13012197 1.0000000
## M85-M31    0.677847496 -0.562554226  1.91824922 0.9226687
## M86-M31   -0.050272482 -1.920248437  1.81970347 1.0000000
## M89-M31   -0.503743967 -1.657640705  0.65015277 0.9925576
## M90-M31   -0.312754015 -1.471535732  0.84602770 0.9999931
## R12-M31   -0.958131261 -2.293561206  0.37729868 0.5397444
## M44-M34    0.294957546 -1.002285626  1.59220072 0.9999996
## M67-M34   -1.528538069 -2.927899938 -0.12917620 0.0164009
## M70-M34   -0.737455621 -2.200943294  0.72603205 0.9638644
## M71-M34   -0.898761207 -2.236203859  0.43868144 0.6653679
## M72-M34    0.046413015 -1.208004990  1.30083102 1.0000000
## M78-M34   -1.114482925 -3.327068313  1.09810246 0.9640084
## M79-M34   -0.216318447 -1.525301615  1.09266472 1.0000000
## M83-M34   -0.394367440 -1.716640118  0.92790524 0.9999652
## M84-M34    0.101481105 -1.235961546  1.43892376 1.0000000
## M85-M34    0.919605133 -0.389378035  2.22858830 0.5814412
## M86-M34    0.191485155 -1.724670000  2.10764031 1.0000000
## M89-M34   -0.261986330 -1.489307879  0.96533522 0.9999999
## M90-M34   -0.070996378 -1.302911786  1.16091903 1.0000000
## R12-M34   -0.716373624 -2.115735493  0.68298824 0.9577603
## M67-M44   -1.823495615 -2.915872124 -0.73111911 0.0000013
## M70-M44   -1.032413167 -2.205813766  0.14098743 0.1697352
## M71-M44   -1.193718754 -2.205561450 -0.18187606 0.0050542
## M72-M44   -0.248544531 -1.147787173  0.65069811 0.9999898
## M78-M44   -1.409440471 -3.441829926  0.62294898 0.6065967
## M79-M44   -0.511275993 -1.485190636  0.46263865 0.9462644
## M83-M44   -0.689324986 -1.681029495  0.30237952 0.6021533
## M84-M44   -0.193476441 -1.205319137  0.81836625 1.0000000
## M85-M44    0.624647587 -0.349267056  1.59856223 0.7443888
## M86-M44   -0.103472392 -1.808383936  1.60143915 1.0000000
## M89-M44   -0.556943876 -1.417984515  0.30409676 0.7308797
## M90-M44   -0.365953924 -1.233530070  0.50162222 0.9951070
## R12-M44   -1.011331170 -2.103707679  0.08104534 0.1099401
## M70-M67    0.791082449 -0.494313506  2.07647840 0.8045262
## M71-M67    0.629776862 -0.510047438  1.76960116 0.9149291
## M72-M67    1.574951085  0.533792537  2.61610963 0.0000238
## M78-M67    0.414055145 -1.684987659  2.51309795 1.0000000
## M79-M67    1.312219623  0.205926928  2.41851232 0.0046177
## M83-M67    1.134170630  0.012185057  2.25615620 0.0441804
## M84-M67    1.630019174  0.490194874  2.76984347 0.0001009
## M85-M67    2.448143203  1.341850508  3.55443590 0.0000000
## M86-M67    1.720023224 -0.063820145  3.50386659 0.0740531
## M89-M67    1.266551739  0.258204192  2.27489929 0.0016601
## M90-M67    1.457541691  0.443607689  2.47147569 0.0000888
## R12-M67    0.812164445 -0.399718483  2.02404737 0.6702812
## M71-M70   -0.161305587 -1.379000936  1.05638976 1.0000000
## M72-M70    0.783868636 -0.342005980  1.90974325 0.5989881
## M78-M70   -0.377027304 -2.519353894  1.76529929 1.0000000
## M79-M70    0.521137174 -0.665229571  1.70750392 0.9919860
## M83-M70    0.343088181 -0.857925622  1.54410198 0.9999825
## M84-M70    0.838936726 -0.378758624  2.05663207 0.6191409
## M85-M70    1.657060754  0.470694009  2.84342750 0.0001770
## M86-M70    0.928940775 -0.905638112  2.76351966 0.9620775
## M89-M70    0.475469291 -0.620134335  1.57107292 0.9930695
## M90-M70    0.666459243 -0.434288090  1.76720658 0.8255764
## R12-M70    0.021081997 -1.264313958  1.30647795 1.0000000
## M72-M71    0.945174223 -0.011147024  1.90149547 0.0569698
## M78-M71   -0.215721717 -2.274002528  1.84255909 1.0000000
## M79-M71    0.682442761 -0.344408119  1.70929364 0.6850927
## M83-M71    0.504393768 -0.539345103  1.54813264 0.9765164
## M84-M71    1.000242312 -0.062649328  2.06313395 0.0942988
## M85-M71    1.818366341  0.791515461  2.84521722 0.0000002
## M86-M71    1.090246362 -0.645448418  2.82594114 0.7760160
## M89-M71    0.636774877 -0.283716285  1.55726604 0.6113393
## M90-M71    0.827764829 -0.098842620  1.75437228 0.1501077
## R12-M71    0.182387583 -0.957436717  1.32221188 1.0000000
## M78-M72   -1.160895940 -3.166221745  0.84442986 0.8752333
## M79-M72   -0.262731462 -1.178828846  0.65336592 0.9999814
## M83-M72   -0.440780455 -1.375768433  0.49420752 0.9820216
## M84-M72    0.055068090 -0.901253157  1.01138934 1.0000000
## M85-M72    0.873192118 -0.042905266  1.78928950 0.0832143
## M86-M72    0.145072139 -1.527485202  1.81762948 1.0000000
## M89-M72   -0.308399346 -1.103456367  0.48665768 0.9983433
## M90-M72   -0.117409393 -0.919539715  0.68472093 1.0000000
## R12-M72   -0.762786640 -1.803945188  0.27837191 0.4976327
## M79-M78    0.898164478 -1.141738458  2.93806741 0.9917606
## M83-M78    0.720115485 -1.328340540  2.76857151 0.9995792
## M84-M78    1.215964030 -0.842316781  3.27424484 0.8542333
## M85-M78    2.034088058 -0.005814878  4.07399099 0.0516329
## M86-M78    1.305968079 -1.167777588  3.77971375 0.9433674
## M89-M78    0.852496595 -1.135991584  2.84098477 0.9940237
## M90-M78    1.043486547 -0.947840303  3.03481340 0.9471748
## R12-M78    0.398109300 -1.700933503  2.49715210 1.0000000
## M83-M79   -0.178048993 -1.185061863  0.82896388 1.0000000
## M84-M79    0.317799552 -0.709051328  1.34465043 0.9999378
## M85-M79    1.135923580  0.146425313  2.12542185 0.0078509
## M86-M79    0.407803601 -1.306057671  2.12166487 0.9999992
## M89-M79   -0.045667883 -0.924296400  0.83296063 1.0000000
## M90-M79    0.145322069 -0.739712087  1.03035622 1.0000000
## R12-M79   -0.500055178 -1.606347872  0.60623752 0.9887970
## M84-M83    0.495848544 -0.547890326  1.53958742 0.9804490
## M85-M83    1.313972573  0.306959702  2.32098544 0.0007844
## M86-M83    0.585852594 -1.138180051  2.30988524 0.9997435
## M89-M83    0.132381109 -0.765926283  1.03068850 1.0000000
## M90-M83    0.323371062 -0.581202626  1.22794475 0.9994648
## R12-M83   -0.322006185 -1.443991758  0.79997939 0.9999812
## M85-M84    0.818124028 -0.208726852  1.84497491 0.3315983
## M86-M84    0.090004050 -1.645690731  1.82569883 1.0000000
## M89-M84   -0.363467435 -1.283958597  0.55702373 0.9979006
## M90-M84   -0.172477483 -1.099084933  0.75412997 1.0000000
## R12-M84   -0.817854729 -1.957679029  0.32196957 0.5395860
## M86-M85   -0.728119979 -2.441981251  0.98574129 0.9946538
## M89-M85   -1.181591463 -2.060219980 -0.30296295 0.0004153
## M90-M85   -0.990601511 -1.875635667 -0.10556736 0.0115438
## R12-M85   -1.635978758 -2.742271452 -0.52968606 0.0000434
## M89-M86   -0.453471485 -2.105803688  1.19886072 0.9999909
## M90-M86   -0.262481533 -1.918228825  1.39326576 1.0000000
## R12-M86   -0.907858779 -2.691702148  0.87598459 0.9601211
## M90-M89    0.190989952 -0.568066505  0.95004641 0.9999979
## R12-M89   -0.454387294 -1.462734841  0.55396025 0.9891869
## R12-M90   -0.645377246 -1.659311249  0.36855676 0.7562254
## 
## $Species
##                diff         lwr       upr     p adj
## OLI-CIR  0.49329645  0.05844870 0.9281442 0.0216821
## VUL-CIR  0.49150464 -0.01501893 0.9980282 0.0593467
## VUL-OLI -0.00179181 -0.31398055 0.3103969 0.9998990
## 
## $ReproductiveModeSimple
##                diff         lwr       upr     p adj
## NR-ASEX  0.07527977 -0.22845386 0.3790134 0.8285485
## SEX-ASEX 0.22784165 -0.02623834 0.4819216 0.0889673
## SEX-NR   0.15256187 -0.18618096 0.4913047 0.5385271
TukeyHSD(aov(chaoMod))
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = chaoMod)
## 
## $PopID
##                   diff          lwr          upr     p adj
## M108-M107  0.479950567 -0.350709648  1.310610783 0.8772257
## M109-M107  0.645872725 -0.204335111  1.496080561 0.4246334
## M26-M107   0.756218971 -0.117286544  1.629724486 0.1923080
## M28-M107   0.386702299 -0.343732247  1.117136846 0.9419521
## M31-M107   0.044191388 -0.892692017  0.981074793 1.0000000
## M34-M107  -0.396043363 -1.377778809  0.585692083 0.9972282
## M44-M107   0.316236045 -0.450130942  1.082603032 0.9962950
## M67-M107  -0.869586527 -1.719794363 -0.019378691 0.0384702
## M70-M107  -0.954338319 -1.856119908 -0.052556729 0.0251345
## M71-M107  -0.805865845 -1.605520275 -0.006211415 0.0457733
## M72-M107   0.476048303 -0.254386243  1.206482849 0.7186610
## M78-M107  -1.763127646 -3.235730816 -0.290524477 0.0039359
## M79-M107   0.147247991 -0.628882026  0.923378008 1.0000000
## M83-M107   0.096387857 -0.690751647  0.883527361 1.0000000
## M84-M107   0.459273512 -0.340380918  1.258927942 0.8829182
## M85-M107   0.977481399  0.201351381  1.753611416 0.0015755
## M86-M107  -0.095063331 -1.346535380  1.156408718 1.0000000
## M89-M107   0.010412917 -0.697002763  0.717828597 1.0000000
## M90-M107  -0.129082098 -0.840417009  0.582252813 1.0000000
## R12-M107  -0.541320764 -1.391528600  0.308887073 0.7559491
## M109-M108  0.165922157 -0.664738058  0.996582373 1.0000000
## M26-M108   0.276268404 -0.578222622  1.130759430 0.9998786
## M28-M108  -0.093248268 -0.800834062  0.614337525 1.0000000
## M31-M108  -0.435759180 -1.354940089  0.483421729 0.9809226
## M34-M108  -0.875993930 -1.840850145  0.088862284 0.1309101
## M44-M108  -0.163714522 -0.908336163  0.580907118 0.9999998
## M67-M108  -1.349537095 -2.180197310 -0.518876879 0.0000030
## M70-M108  -1.434288886 -2.317664834 -0.550912939 0.0000031
## M71-M108  -1.285816413 -2.064655443 -0.506977382 0.0000018
## M72-M108  -0.003902265 -0.711488058  0.703683529 1.0000000
## M78-M108  -2.243078214 -3.704482717 -0.781673710 0.0000159
## M79-M108  -0.332702577 -1.087368617  0.421963464 0.9916513
## M83-M108  -0.383562711 -1.149546801  0.382421380 0.9661265
## M84-M108  -0.020677055 -0.799516086  0.758161976 1.0000000
## M85-M108   0.497530831 -0.257135209  1.252196872 0.6992385
## M86-M108  -0.575013899 -1.813289037  0.663261240 0.9848130
## M89-M108  -0.469537650 -1.153335847  0.214260546 0.6256382
## M90-M108  -0.609032665 -1.296884673  0.078819342 0.1614578
## R12-M108  -1.021271331 -1.851931546 -0.190611116 0.0024650
## M26-M109   0.110346247 -0.763159269  0.983851762 1.0000000
## M28-M109  -0.259170425 -0.989604972  0.471264121 0.9995205
## M31-M109  -0.601681337 -1.538564742  0.335202068 0.7424993
## M34-M109  -1.041916088 -2.023651534 -0.060180642 0.0242174
## M44-M109  -0.329636680 -1.096003667  0.436730307 0.9937919
## M67-M109  -1.515459252 -2.365667088 -0.665251416 0.0000001
## M70-M109  -1.600211044 -2.501992633 -0.698429454 0.0000002
## M71-M109  -1.451738570 -2.251393000 -0.652084140 0.0000001
## M72-M109  -0.169824422 -0.900258968  0.560610124 0.9999995
## M78-M109  -2.409000371 -3.881603540 -0.936397202 0.0000024
## M79-M109  -0.498624734 -1.274754751  0.277505284 0.7419226
## M83-M109  -0.549484868 -1.336624372  0.237654636 0.5940273
## M84-M109  -0.186599212 -0.986253642  0.613055217 0.9999994
## M85-M109   0.331608674 -0.444521343  1.107738691 0.9942796
## M86-M109  -0.740936056 -1.992408105  0.510535993 0.8519785
## M89-M109  -0.635459808 -1.342875488  0.071955873 0.1433780
## M90-M109  -0.774954823 -1.486289734 -0.063619912 0.0170089
## R12-M109  -1.187193488 -2.037401324 -0.336985652 0.0001777
## M28-M26   -0.369516672 -1.126942012  0.387908669 0.9740919
## M31-M26   -0.712027584 -1.670103277  0.246048110 0.4688492
## M34-M26   -1.152262334 -2.154241878 -0.150282790 0.0076334
## M44-M26   -0.439982926 -1.232117308  0.352151456 0.9111155
## M67-M26   -1.625805499 -2.499311014 -0.752299983 0.0000000
## M70-M26   -1.710557290 -2.634336785 -0.786777795 0.0000000
## M71-M26   -1.562084816 -2.386466844 -0.737702789 0.0000000
## M72-M26   -0.280170668 -1.037596009  0.477254672 0.9991365
## M78-M26   -2.519346618 -4.005522453 -1.033170782 0.0000008
## M79-M26   -0.608970980 -1.410554615  0.192612655 0.4245138
## M83-M26   -0.659831114 -1.472079302  0.152417073 0.2960066
## M84-M26   -0.296945459 -1.121327487  0.527436569 0.9994055
## M85-M26    0.221262427 -0.580321208  1.022846063 0.9999900
## M86-M26   -0.851282302 -2.118697313  0.416132709 0.6664356
## M89-M26   -0.745806054 -1.481058028 -0.010554080 0.0424354
## M90-M26   -0.885301069 -1.624324665 -0.146277473 0.0038995
## R12-M26   -1.297539735 -2.171045250 -0.424034220 0.0000385
## M31-M28   -0.342510912 -1.172228801  0.487206978 0.9962766
## M34-M28   -0.782745662 -1.662794381  0.097303057 0.1556335
## M44-M28   -0.070466254 -0.701338365  0.560405857 1.0000000
## M67-M28   -1.256288827 -1.986723373 -0.525854280 0.0000005
## M70-M28   -1.341040618 -2.130908517 -0.551172720 0.0000007
## M71-M28   -1.192568144 -1.863484286 -0.521652003 0.0000002
## M72-M28    0.089346004 -0.497353142  0.676045150 1.0000000
## M78-M28   -2.149829946 -3.556685075 -0.742974816 0.0000180
## M79-M28   -0.239454308 -0.882151022  0.403242405 0.9990456
## M83-M28   -0.290314443 -0.946264029  0.365635144 0.9912561
## M84-M28    0.072571213 -0.598344929  0.743487354 1.0000000
## M85-M28    0.590779099 -0.051917614  1.233475813 0.1173364
## M86-M28   -0.481765630 -1.655163922  0.691632662 0.9965246
## M89-M28   -0.376289382 -0.934069094  0.181490329 0.6584873
## M90-M28   -0.515784397 -1.078526449  0.046957654 0.1204262
## R12-M28   -0.928023063 -1.658457609 -0.197588517 0.0013314
## M34-M31   -0.440234751 -1.497921678  0.617452176 0.9958675
## M44-M31    0.272044657 -0.589474645  1.133563960 0.9999154
## M67-M31   -0.913777915 -1.850661320  0.023105490 0.0657153
## M70-M31   -0.998529707 -1.982454086 -0.014605327 0.0421960
## M71-M31   -0.850057233 -1.741317216  0.041202750 0.0827010
## M72-M31    0.431856915 -0.397860974  1.261574805 0.9505133
## M78-M31   -1.807319034 -3.331608129 -0.283029940 0.0046437
## M79-M31    0.103056603 -0.767158861  0.973272067 1.0000000
## M83-M31    0.052196469 -0.827852250  0.932245188 1.0000000
## M84-M31    0.415082125 -0.476177858  1.306342107 0.9843104
## M85-M31    0.933290011  0.063074547  1.803505475 0.0210765
## M86-M31   -0.139254719 -1.451153892  1.172644454 1.0000000
## M89-M31   -0.033778471 -0.843305551  0.775748610 1.0000000
## M90-M31   -0.173273486 -0.986227669  0.639680698 0.9999999
## R12-M31   -0.585512151 -1.522395556  0.351371254 0.7834814
## M44-M34    0.712279408 -0.197813708  1.622372524 0.3657738
## M67-M34   -0.473543164 -1.455278610  0.508192282 0.9770166
## M70-M34   -0.558294956 -1.585018461  0.468428550 0.9260670
## M71-M34   -0.409822482 -1.348117920  0.528472956 0.9925289
## M72-M34    0.872091666 -0.007957053  1.752140385 0.0553397
## M78-M34   -1.367084283 -2.919344318  0.185175751 0.1684780
## M79-M34    0.543291354 -0.375038067  1.461620775 0.8528054
## M83-M34    0.492431220 -0.435221582  1.420084021 0.9404971
## M84-M34    0.855316875 -0.082978563  1.793612313 0.1263885
## M85-M34    1.373524762  0.455195341  2.291854183 0.0000321
## M86-M34    0.300980032 -1.043316592  1.645276655 0.9999997
## M89-M34    0.406456280 -0.454582666  1.267495226 0.9817910
## M90-M34    0.266961265 -0.597300545  1.131223075 0.9999398
## R12-M34   -0.145277401 -1.127012847  0.836458046 1.0000000
## M67-M44   -1.185822572 -1.952189559 -0.419455585 0.0000127
## M70-M44   -1.270574364 -2.093784562 -0.447364166 0.0000136
## M71-M44   -1.122101890 -1.831969625 -0.412234155 0.0000069
## M72-M44    0.159812258 -0.471059853  0.790684369 0.9999976
## M78-M44   -2.079363692 -3.505205579 -0.653521804 0.0000618
## M79-M44   -0.168988054 -0.852247008  0.514270900 0.9999984
## M83-M44   -0.219848188 -0.915587789  0.475891413 0.9999145
## M84-M44    0.143037467 -0.566830268  0.852905202 1.0000000
## M85-M44    0.661245353 -0.022013600  1.344504307 0.0712976
## M86-M44   -0.411299376 -1.607396063  0.784797311 0.9996955
## M89-M44   -0.305823128 -0.909894265  0.298248009 0.9621934
## M90-M44   -0.445318143 -1.053974326  0.163338040 0.5005306
## R12-M44   -0.857556809 -1.623923796 -0.091189822 0.0115822
## M70-M67   -0.084751792 -0.986533381  0.817029798 1.0000000
## M71-M67    0.063720682 -0.735933748  0.863375112 1.0000000
## M72-M67    1.345634830  0.615200284  2.076069376 0.0000000
## M78-M67   -0.893541119 -2.366144288  0.579062050 0.8230199
## M79-M67    1.016834518  0.240704501  1.792964536 0.0007203
## M83-M67    0.965974384  0.178834880  1.753113888 0.0025490
## M84-M67    1.328860040  0.529205610  2.128514469 0.0000015
## M85-M67    1.847067926  1.070937909  2.623197943 0.0000000
## M86-M67    0.774523196 -0.476948853  2.025995245 0.7971179
## M89-M67    0.879999444  0.172583764  1.587415125 0.0019879
## M90-M67    0.740504429  0.029169519  1.451839340 0.0310015
## R12-M67    0.328265764 -0.521942072  1.178473600 0.9984472
## M71-M70    0.148472474 -0.705813121  1.002758069 1.0000000
## M72-M70    1.430386622  0.640518723  2.220254520 0.0000001
## M78-M70   -0.808789328 -2.311758643  0.694179988 0.9326269
## M79-M70    1.101586310  0.269279591  1.933893029 0.0005827
## M83-M70    1.050726176  0.208143676  1.893308675 0.0018988
## M84-M70    1.413611831  0.559326236  2.267897426 0.0000017
## M85-M70    1.931819717  1.099512999  2.764126436 0.0000000
## M86-M70    0.859274988 -0.427791040  2.146341016 0.6772477
## M89-M70    0.964751236  0.196120235  1.733382237 0.0016806
## M90-M70    0.825256221  0.053016604  1.597495838 0.0221082
## R12-M70    0.413017555 -0.488764034  1.314799145 0.9869993
## M72-M71    1.281914148  0.610998007  1.952830289 0.0000000
## M78-M71   -0.957261801 -2.401268012  0.486744409 0.6896892
## M79-M71    0.953113836  0.232716969  1.673510703 0.0005873
## M83-M71    0.902253702  0.170008906  1.634498498 0.0023685
## M84-M71    1.265139357  0.519457757  2.010820958 0.0000008
## M85-M71    1.783347244  1.062950376  2.503744111 0.0000000
## M86-M71    0.710802514 -0.506890441  1.928495469 0.8670571
## M89-M71    0.816278762  0.170499553  1.462057972 0.0014696
## M90-M71    0.676783747  0.026713599  1.326853895 0.0309693
## R12-M71    0.264545082 -0.535109348  1.064199511 0.9998283
## M78-M72   -2.239175949 -3.646031079 -0.832320820 0.0000056
## M79-M72   -0.328800312 -0.971497025  0.313896401 0.9580889
## M83-M72   -0.379660446 -1.035610033  0.276289140 0.8755448
## M84-M72   -0.016774791 -0.687690932  0.654141351 1.0000000
## M85-M72    0.501433096 -0.141263618  1.144129809 0.3718473
## M86-M72   -0.571111634 -1.744509926  0.602286658 0.9747216
## M89-M72   -0.465635386 -1.023415097  0.092144326 0.2488381
## M90-M72   -0.605130401 -1.167872453 -0.042388349 0.0203386
## R12-M72   -1.017369066 -1.747803613 -0.286934520 0.0001885
## M79-M78    1.910375637  0.479262597  3.341488677 0.0004861
## M83-M78    1.859515503  0.422401963  3.296629043 0.0009282
## M84-M78    2.222401159  0.778394948  3.666407369 0.0000148
## M85-M78    2.740609045  1.309496005  4.171722085 0.0000000
## M86-M78    1.668064315 -0.067415163  3.403543794 0.0765155
## M89-M78    1.773540564  0.378498029  3.168583098 0.0013151
## M90-M78    1.634045548  0.237011517  3.031079580 0.0058174
## R12-M78    1.221806883 -0.250796286  2.694410052 0.2591282
## M83-M79   -0.050860134 -0.757339460  0.655619191 1.0000000
## M84-M79    0.312025521 -0.408371346  1.032422389 0.9932454
## M85-M79    0.830233408  0.136041616  1.524425199 0.0040107
## M86-M79   -0.242311322 -1.444686775  0.960064131 1.0000000
## M89-M79   -0.136835074 -0.753245152  0.479575004 0.9999998
## M90-M79   -0.276330089 -0.897234103  0.344573925 0.9906668
## R12-M79   -0.688568754 -1.464698772  0.087561263 0.1588752
## M84-M83    0.362885655 -0.369359140  1.095130451 0.9695160
## M85-M83    0.881093542  0.174614216  1.587572867 0.0018951
## M86-M83   -0.191451188 -1.400962463  1.018060087 1.0000000
## M89-M83   -0.085974940 -0.716190917  0.544241038 1.0000000
## M90-M83   -0.225469955 -0.860082111  0.409142201 0.9995113
## R12-M83   -0.637708620 -1.424848124  0.149430883 0.3008788
## M85-M84    0.518207886 -0.202188981  1.238604754 0.5347069
## M86-M84   -0.554336843 -1.772029798  0.663356112 0.9878667
## M89-M84   -0.448860595 -1.094639805  0.196918614 0.6023924
## M90-M84   -0.588355610 -1.238425759  0.061714538 0.1344962
## R12-M84   -1.000594276 -1.800248706 -0.200939846 0.0017815
## M86-M85   -1.072544730 -2.274920182  0.129830723 0.1519716
## M89-M85   -0.967068482 -1.583478559 -0.350658404 0.0000086
## M90-M85   -1.106563497 -1.727467511 -0.485659483 0.0000001
## R12-M85   -1.518802162 -2.294932180 -0.742672145 0.0000000
## M89-M86    0.105476248 -1.053732909  1.264685405 1.0000000
## M90-M86   -0.034018767 -1.195623812  1.127586278 1.0000000
## R12-M86   -0.446257432 -1.697729482  0.805214617 0.9994851
## M90-M89   -0.139495015 -0.672018193  0.393028163 0.9999958
## R12-M89   -0.551733681 -1.259149361  0.155682000 0.3725295
## R12-M90   -0.412238666 -1.123573576  0.299096245 0.8743051
## 
## $Species
##               diff         lwr        upr     p adj
## OLI-CIR  0.3610004  0.05591592 0.66608496 0.0156291
## VUL-CIR  0.2098717 -0.14549980 0.56524321 0.3463080
## VUL-OLI -0.1511287 -0.37015701 0.06789953 0.2361629
# do some nice plots
models <- list("Simpson" = simpMod ,"Shannon" = shanMod,"Chao1" = chaoMod)
effecSizePlots <- list()
for(vrb in c("Species","PopID")){
    for(n in names(models)){
    mod  <- models[[n]]
    prd  <- make_predictions(mod, pred = vrb, data = alphaDiv)
    effecSizePlots[[paste(n,vrb,sep="_")]] <-
        ggplot(prd, aes_string(y = vrb, x=n ,xmin = "ymin", xmax = "ymax")) +
        geom_point() +
        geom_errorbar() +
        ggtitle(paste0("Effect size of ",n))
    }
}

pp <- plot_grid(plotlist=effecSizePlots, ncol = 3, rel_heights =c(0.3,0.7))
print(pp)

This last part is has to be taken carefully, as the design is unbalanced and Anova calculates Type II anovas, so only main effects are considered, which is fine as there are no interactions.